mesh_bsurfaces.py

                    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
        me = object.data
        bm = bmesh.new()
        bm.from_mesh(me)

        num_faces_created = 0
        for i in range(len(faces_verts_idx)):
            if i not in faces_to_not_include_idx:
                bm.faces.new([bm.verts[v] for v in faces_verts_idx[i]])

                num_faces_created += 1

        bm.to_mesh(me)
        bm.free()

        for v_idx in selected_verts_idx:
            self.main_object.data.vertices[v_idx].select = True

        bpy.ops.object.mode_set('INVOKE_REGION_WIN', mode='EDIT')
        bpy.ops.mesh.normals_make_consistent(inside=False)
        bpy.ops.object.mode_set('INVOKE_REGION_WIN', mode='OBJECT')

        self.update()

        return num_faces_created

    # Crosshatch skinning
    def crosshatch_surface_invoke(self, ob_original_splines):
        self.is_crosshatch = False
        self.crosshatch_merge_distance = 0

        objects_to_delete = []  # duplicated strokes to be deleted.

        # If the main object uses modifiers deactivate them temporarily until the surface is joined
        # (without this the surface verts merging with the main object doesn't work well)
        self.modifiers_prev_viewport_state = []
        if len(self.main_object.modifiers) > 0:
            for m_idx in range(len(self.main_object.modifiers)):
                self.modifiers_prev_viewport_state.append(
                                    self.main_object.modifiers[m_idx].show_viewport
                                    )
                self.main_object.modifiers[m_idx].show_viewport = False

        bpy.ops.object.select_all('INVOKE_REGION_WIN', action='DESELECT')
        ob_original_splines.select_set(True)
        bpy.context.view_layer.objects.active = ob_original_splines

        if len(ob_original_splines.data.splines) >= 2:
            bpy.ops.object.duplicate('INVOKE_REGION_WIN')
            ob_splines = bpy.context.object
            ob_splines.name = "SURFSKIO_NE_STR"

            # Get estimative merge distance (sum up the distances from the first point to
            # all other points, then average them and then divide them)
            first_point_dist_sum = 0
            first_dist = 0
            second_dist = 0
            coords_first_pt = ob_splines.data.splines[0].bezier_points[0].co
            for i in range(len(ob_splines.data.splines)):
                sp = ob_splines.data.splines[i]

                if coords_first_pt != sp.bezier_points[0].co:
                    first_dist = (coords_first_pt - sp.bezier_points[0].co).length

                if coords_first_pt != sp.bezier_points[len(sp.bezier_points) - 1].co:
                    second_dist = (coords_first_pt - sp.bezier_points[len(sp.bezier_points) - 1].co).length

                first_point_dist_sum += first_dist + second_dist

                if i == 0:
                    if first_dist != 0:
                        shortest_dist = first_dist
                    elif second_dist != 0:
                        shortest_dist = second_dist

                if shortest_dist > first_dist and first_dist != 0:
                    shortest_dist = first_dist

                if shortest_dist > second_dist and second_dist != 0:
                    shortest_dist = second_dist

            self.crosshatch_merge_distance = shortest_dist / 20

            # Recalculation of merge distance

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

            ob_calc_merge_dist = bpy.context.object
            ob_calc_merge_dist.name = "SURFSKIO_CALC_TMP"

            objects_to_delete.append(ob_calc_merge_dist)

            # Smooth out strokes a little to improve crosshatch detection
            bpy.ops.object.editmode_toggle('INVOKE_REGION_WIN')
            bpy.ops.curve.select_all('INVOKE_REGION_WIN', action='SELECT')

            for i in range(4):
                bpy.ops.curve.smooth('INVOKE_REGION_WIN')

            bpy.ops.curve.select_all('INVOKE_REGION_WIN', action='DESELECT')
            bpy.ops.object.editmode_toggle('INVOKE_REGION_WIN')

            # Convert curves into mesh
            ob_calc_merge_dist.data.resolution_u = 12
            bpy.ops.object.convert(target='MESH', keep_original=False)

            # Find "intersection-nodes"
            bpy.ops.object.editmode_toggle('INVOKE_REGION_WIN')
            bpy.ops.mesh.select_all('INVOKE_REGION_WIN', action='SELECT')
            bpy.ops.mesh.remove_doubles('INVOKE_REGION_WIN',
                                        threshold=self.crosshatch_merge_distance)
            bpy.ops.mesh.select_all('INVOKE_REGION_WIN', action='DESELECT')
            bpy.ops.object.editmode_toggle('INVOKE_REGION_WIN')

            # Remove verts with less than three edges
            verts_edges_count = {}
            for ed in ob_calc_merge_dist.data.edges:
                v = ed.vertices

                if v[0] not in verts_edges_count:
                    verts_edges_count[v[0]] = 0

                if v[1] not in verts_edges_count:
                    verts_edges_count[v[1]] = 0

                verts_edges_count[v[0]] += 1
                verts_edges_count[v[1]] += 1

            nodes_verts_coords = []
            for v_idx in verts_edges_count:
                v = ob_calc_merge_dist.data.vertices[v_idx]

                if verts_edges_count[v_idx] < 3:
                    v.select = True

            # Remove them
            bpy.ops.object.editmode_toggle('INVOKE_REGION_WIN')
            bpy.ops.mesh.delete('INVOKE_REGION_WIN', type='VERT')
            bpy.ops.mesh.select_all('INVOKE_REGION_WIN', action='SELECT')

            # Remove doubles to discard very near verts from calculations of distance
            bpy.ops.mesh.remove_doubles(
                        'INVOKE_REGION_WIN',
                        threshold=self.crosshatch_merge_distance * 4.0
                        )
            bpy.ops.object.editmode_toggle('INVOKE_REGION_WIN')

            # Get all coords of the resulting nodes
            nodes_verts_coords = [(v.co[0], v.co[1], v.co[2]) for
                                   v in ob_calc_merge_dist.data.vertices]

            # Check if the strokes are a crosshatch
            if len(nodes_verts_coords) >= 3:
                self.is_crosshatch = True

                shortest_dist = None
                for co_1 in nodes_verts_coords:
                    for co_2 in nodes_verts_coords:
                        if co_1 != co_2:
                            dist = (Vector(co_1) - Vector(co_2)).length

                            if shortest_dist is not None:
                                if dist < shortest_dist:
                                    shortest_dist = dist
                            else:
                                shortest_dist = dist

                self.crosshatch_merge_distance = shortest_dist / 3

            bpy.ops.object.select_all('INVOKE_REGION_WIN', action='DESELECT')
            ob_splines.select_set(True)
            bpy.context.view_layer.objects.active = ob_splines

            # Deselect all points
            bpy.ops.object.editmode_toggle('INVOKE_REGION_WIN')
            bpy.ops.curve.select_all('INVOKE_REGION_WIN', action='DESELECT')
            bpy.ops.object.editmode_toggle('INVOKE_REGION_WIN')

            # Smooth splines in a localized way, to eliminate "saw-teeth"
            # like shapes when there are many points
            for sp in ob_splines.data.splines:
                angle_sum = 0

                angle_limit = 2  # Degrees
                for t in range(len(sp.bezier_points)):
                    # Because on each iteration it checks the "next two points"
                    # of the actual. This way it doesn't go out of range
                    if t <= len(sp.bezier_points) - 3:
                        p1 = sp.bezier_points[t]
                        p2 = sp.bezier_points[t + 1]
                        p3 = sp.bezier_points[t + 2]

                        vec_1 = p1.co - p2.co
                        vec_2 = p2.co - p3.co

                        if p2.co != p1.co and p2.co != p3.co:
                            angle = vec_1.angle(vec_2)
                            angle_sum += degrees(angle)

                            if angle_sum >= angle_limit:  # If sum of angles is grater than the limit
                                if (p1.co - p2.co).length <= self.crosshatch_merge_distance:
                                    p1.select_control_point = True
                                    p1.select_left_handle = True
                                    p1.select_right_handle = True

                                    p2.select_control_point = True
                                    p2.select_left_handle = True
                                    p2.select_right_handle = True

                                if (p1.co - p2.co).length <= self.crosshatch_merge_distance:
                                    p3.select_control_point = True
                                    p3.select_left_handle = True
                                    p3.select_right_handle = True

                                angle_sum = 0

                sp.bezier_points[0].select_control_point = False
                sp.bezier_points[0].select_left_handle = False
                sp.bezier_points[0].select_right_handle = False

                sp.bezier_points[len(sp.bezier_points) - 1].select_control_point = False
                sp.bezier_points[len(sp.bezier_points) - 1].select_left_handle = False
                sp.bezier_points[len(sp.bezier_points) - 1].select_right_handle = False

            # Smooth out strokes a little to improve crosshatch detection
            bpy.ops.object.editmode_toggle('INVOKE_REGION_WIN')

            for i in range(15):
                bpy.ops.curve.smooth('INVOKE_REGION_WIN')

            bpy.ops.curve.select_all('INVOKE_REGION_WIN', action='DESELECT')
            bpy.ops.object.editmode_toggle('INVOKE_REGION_WIN')

            # Simplify the splines
            for sp in ob_splines.data.splines:
                angle_sum = 0

                sp.bezier_points[0].select_control_point = True
                sp.bezier_points[0].select_left_handle = True
                sp.bezier_points[0].select_right_handle = True

                sp.bezier_points[len(sp.bezier_points) - 1].select_control_point = True
                sp.bezier_points[len(sp.bezier_points) - 1].select_left_handle = True
                sp.bezier_points[len(sp.bezier_points) - 1].select_right_handle = True

                angle_limit = 15  # Degrees
                for t in range(len(sp.bezier_points)):
                    # Because on each iteration it checks the "next two points"
                    # of the actual. This way it doesn't go out of range
                    if t <= len(sp.bezier_points) - 3:
                        p1 = sp.bezier_points[t]
                        p2 = sp.bezier_points[t + 1]
                        p3 = sp.bezier_points[t + 2]

                        vec_1 = p1.co - p2.co
                        vec_2 = p2.co - p3.co

                        if p2.co != p1.co and p2.co != p3.co:
                            angle = vec_1.angle(vec_2)
                            angle_sum += degrees(angle)
                            # If sum of angles is grater than the limit
                            if angle_sum >= angle_limit:
                                p1.select_control_point = True
                                p1.select_left_handle = True
                                p1.select_right_handle = True

                                p2.select_control_point = True
                                p2.select_left_handle = True
                                p2.select_right_handle = True

                                p3.select_control_point = True
                                p3.select_left_handle = True
                                p3.select_right_handle = True

                                angle_sum = 0

            bpy.ops.object.editmode_toggle('INVOKE_REGION_WIN')
            bpy.ops.curve.select_all(action='INVERT')

            bpy.ops.curve.delete(type='VERT')
            bpy.ops.object.editmode_toggle('INVOKE_REGION_WIN')

            objects_to_delete.append(ob_splines)

            bpy.ops.object.editmode_toggle('INVOKE_REGION_WIN')
            bpy.ops.curve.select_all('INVOKE_REGION_WIN', action='DESELECT')
            bpy.ops.object.editmode_toggle('INVOKE_REGION_WIN')

            # Check if the strokes are a crosshatch
            if self.is_crosshatch:
                all_points_coords = []
                for i in range(len(ob_splines.data.splines)):
                    all_points_coords.append([])

                    all_points_coords[i] = [Vector((x, y, z)) for
                                            x, y, z in [bp.co for
                                            bp in ob_splines.data.splines[i].bezier_points]]

                all_intersections = []
                checked_splines = []
                for i in range(len(all_points_coords)):

                    for t in range(len(all_points_coords[i]) - 1):
                        bp1_co = all_points_coords[i][t]
                        bp2_co = all_points_coords[i][t + 1]

                        for i2 in range(len(all_points_coords)):
                            if i != i2 and i2 not in checked_splines:
                                for t2 in range(len(all_points_coords[i2]) - 1):
                                    bp3_co = all_points_coords[i2][t2]
                                    bp4_co = all_points_coords[i2][t2 + 1]

                                    intersec_coords = intersect_line_line(
                                                            bp1_co, bp2_co, bp3_co, bp4_co
                                                            )
                                    if intersec_coords is not None:
                                        dist = (intersec_coords[0] - intersec_coords[1]).length

                                        if dist <= self.crosshatch_merge_distance * 1.5:
                                            _temp_co, percent1 = intersect_point_line(
                                                                    intersec_coords[0], bp1_co, bp2_co
                                                                    )
                                            if (percent1 >= -0.02 and percent1 <= 1.02):
                                                _temp_co, percent2 = intersect_point_line(
                                                                        intersec_coords[1], bp3_co, bp4_co
                                                                        )
                                                if (percent2 >= -0.02 and percent2 <= 1.02):
                                                    # Format: spline index, first point index from
                                                    # corresponding segment, percentage from first point of
                                                    # actual segment, coords of intersection point
                                                    all_intersections.append(
                                                            (i, t, percent1,
                                                            ob_splines.matrix_world @ intersec_coords[0])
                                                            )
                                                    all_intersections.append(
                                                            (i2, t2, percent2,
                                                            ob_splines.matrix_world @ intersec_coords[1])
                                                            )

                        checked_splines.append(i)
                # Sort list by spline, then by corresponding first point index of segment,
                # and then by percentage from first point of segment: elements 0 and 1 respectively
                all_intersections.sort(key=operator.itemgetter(0, 1, 2))

                self.crosshatch_strokes_coords = {}
                for i in range(len(all_intersections)):
                    if not all_intersections[i][0] in self.crosshatch_strokes_coords:
                        self.crosshatch_strokes_coords[all_intersections[i][0]] = []

                    self.crosshatch_strokes_coords[all_intersections[i][0]].append(
                                                                        all_intersections[i][3]
                                                                        )  # Save intersection coords
            else:
                self.is_crosshatch = False

        # Delete all duplicates
        bpy.ops.object.delete({"selected_objects": objects_to_delete})

        # If the main object has modifiers, turn their "viewport view status" to
        # what it was before the forced deactivation above
        if len(self.main_object.modifiers) > 0:
            for m_idx in range(len(self.main_object.modifiers)):
                self.main_object.modifiers[m_idx].show_viewport = self.modifiers_prev_viewport_state[m_idx]

        self.update()

        return

    # Part of the Crosshatch process that is repeated when the operator is tweaked
    def crosshatch_surface_execute(self, context):
        # If the main object uses modifiers deactivate them temporarily until the surface is joined
        # (without this the surface verts merging with the main object doesn't work well)
        self.modifiers_prev_viewport_state = []
        if len(self.main_object.modifiers) > 0:
            for m_idx in range(len(self.main_object.modifiers)):
                self.modifiers_prev_viewport_state.append(self.main_object.modifiers[m_idx].show_viewport)

                self.main_object.modifiers[m_idx].show_viewport = False

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

        me_name = "SURFSKIO_STK_TMP"
        me = bpy.data.meshes.new(me_name)

        all_verts_coords = []
        all_edges = []
        for st_idx in self.crosshatch_strokes_coords:
            for co_idx in range(len(self.crosshatch_strokes_coords[st_idx])):
                coords = self.crosshatch_strokes_coords[st_idx][co_idx]

                all_verts_coords.append(coords)

                if co_idx > 0:
                    all_edges.append((len(all_verts_coords) - 2, len(all_verts_coords) - 1))

        me.from_pydata(all_verts_coords, all_edges, [])
        ob = object_utils.object_data_add(context, me)

        bpy.ops.object.select_all('INVOKE_REGION_WIN', action='DESELECT')
        ob.select_set(True)
        bpy.context.view_layer.objects.active = ob

        # Get together each vert and its nearest, to the middle position
        verts = ob.data.vertices
        checked_verts = []
        for i in range(len(verts)):
            shortest_dist = None

            if i not in checked_verts:
                for t in range(len(verts)):
                    if i != t and t not in checked_verts:
                        dist = (verts[i].co - verts[t].co).length

                        if shortest_dist is not None:
                            if dist < shortest_dist:
                                shortest_dist = dist
                                nearest_vert = t
                        else:
                            shortest_dist = dist
                            nearest_vert = t

                middle_location = (verts[i].co + verts[nearest_vert].co) / 2

                verts[i].co = middle_location
                verts[nearest_vert].co = middle_location

                checked_verts.append(i)
                checked_verts.append(nearest_vert)

        # Calculate average length between all the generated edges
        ob = bpy.context.object
        lengths_sum = 0
        for ed in ob.data.edges:
            v1 = ob.data.vertices[ed.vertices[0]]
            v2 = ob.data.vertices[ed.vertices[1]]

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

        edges_count = len(ob.data.edges)
        # possible division by zero here
        average_edge_length = lengths_sum / edges_count if edges_count != 0 else 0.0001

        bpy.ops.object.editmode_toggle('INVOKE_REGION_WIN')
        bpy.ops.mesh.select_all('INVOKE_REGION_WIN', action='SELECT')
        bpy.ops.mesh.remove_doubles('INVOKE_REGION_WIN',
                                    threshold=average_edge_length / 15.0)
        bpy.ops.object.editmode_toggle('INVOKE_REGION_WIN')

        final_points_ob = bpy.context.view_layer.objects.active

        # Make a dictionary with the verts related to each vert
        related_key_verts = {}
        for ed in final_points_ob.data.edges:
            if not ed.vertices[0] in related_key_verts:
                related_key_verts[ed.vertices[0]] = []

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

            if not ed.vertices[1] in related_key_verts[ed.vertices[0]]:
                related_key_verts[ed.vertices[0]].append(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
                        for f_verts in faces_verts_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.
                        for f_verts in faces_verts_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
        for i in range(len(faces_verts_idx)):
            for t in range(len(faces_verts_idx)):
                if i != t:
                    verts_in_common = 0

                    if len(faces_verts_idx[i]) == 4 and len(faces_verts_idx[t]) == 3:
                        for v_idx in faces_verts_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 surface
        all_surface_verts_co = []
        for i in range(len(final_points_ob.data.vertices)):
            coords = final_points_ob.data.vertices[i].co
            all_surface_verts_co.append([coords[0], coords[1], coords[2]])

        # Verts of each face.
        all_surface_faces = []
        for i in range(len(faces_verts_idx)):
            if i not in faces_to_not_include_idx:
                face = []
                for v_idx in faces_verts_idx[i]:
                    face.append(v_idx)

                all_surface_faces.append(face)

        # Build the mesh
        surf_me_name = "SURFSKIO_surface"
        me_surf = bpy.data.meshes.new(surf_me_name)
        me_surf.from_pydata(all_surface_verts_co, [], all_surface_faces)
        ob_surface = object_utils.object_data_add(context, me_surf)

        # Delete final points temporal object
        bpy.ops.object.delete({"selected_objects": [final_points_ob]})

        # Delete isolated verts if there are any
        bpy.ops.object.select_all('INVOKE_REGION_WIN', action='DESELECT')
        ob_surface.select_set(True)
        bpy.context.view_layer.objects.active = ob_surface

        bpy.ops.object.editmode_toggle('INVOKE_REGION_WIN')
        bpy.ops.mesh.select_all(action='DESELECT')
        bpy.ops.mesh.select_face_by_sides(type='NOTEQUAL')
        bpy.ops.mesh.delete()
        bpy.ops.object.editmode_toggle('INVOKE_REGION_WIN')

        # Join crosshatch results with original mesh

        # Calculate a distance to merge the verts of the crosshatch surface to the main object
        edges_length_sum = 0
        for ed in ob_surface.data.edges:
            edges_length_sum += (
                                ob_surface.data.vertices[ed.vertices[0]].co -
                                ob_surface.data.vertices[ed.vertices[1]].co
                                ).length

        # Make dictionary with all the verts connected to each vert, on the new surface object.
        surface_connected_verts = {}
        for ed in ob_surface.data.edges:
            if not ed.vertices[0] in surface_connected_verts:
                surface_connected_verts[ed.vertices[0]] = []

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

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

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

        # Duplicate the new surface object, and use shrinkwrap to
        # calculate later the nearest verts to the main object
        bpy.ops.object.editmode_toggle('INVOKE_REGION_WIN')
        bpy.ops.mesh.select_all('INVOKE_REGION_WIN', action='DESELECT')
        bpy.ops.object.editmode_toggle('INVOKE_REGION_WIN')

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

        final_ob_duplicate = bpy.context.view_layer.objects.active

        bpy.ops.object.modifier_add('INVOKE_REGION_WIN', type='SHRINKWRAP')
        shrinkwrap_modifier = final_ob_duplicate.modifiers[-1]
        shrinkwrap_modifier.wrap_method = "NEAREST_VERTEX"
        shrinkwrap_modifier.target = self.main_object

        bpy.ops.object.modifier_apply('INVOKE_REGION_WIN', modifier=shrinkwrap_modifier.name)

        # Make list with verts of original mesh as index and coords as value
        main_object_verts_coords = []
        for v in self.main_object.data.vertices:
            coords = self.main_object.matrix_world @ v.co

            # To avoid problems when taking "-0.00" as a different value as "0.00"
            for c in range(len(coords)):
                if "%.3f" % coords[c] == "-0.00":
                    coords[c] = 0

            main_object_verts_coords.append(["%.3f" % coords[0], "%.3f" % coords[1], "%.3f" % coords[2]])

        tuple(main_object_verts_coords)

        # Determine which verts will be merged, snap them to the nearest verts
        # on the original verts, and get them selected
        crosshatch_verts_to_merge = []
        if self.automatic_join:
            for i in range(len(ob_surface.data.vertices)-1):
                # Calculate the distance from each of the connected verts to the actual vert,
                # and compare it with the distance they would have if joined.
                # If they don't change much, that vert can be joined
                merge_actual_vert = True
                try:
                    if len(surface_connected_verts[i]) < 4:
                        for c_v_idx in surface_connected_verts[i]:
                            points_original = []
                            points_original.append(ob_surface.data.vertices[c_v_idx].co)
                            points_original.append(ob_surface.data.vertices[i].co)

                            points_target = []
                            points_target.append(ob_surface.data.vertices[c_v_idx].co)
                            points_target.append(final_ob_duplicate.data.vertices[i].co)

                            vec_A = points_original[0] - points_original[1]
                            vec_B = points_target[0] - points_target[1]

                            dist_A = (points_original[0] - points_original[1]).length
                            dist_B = (points_target[0] - points_target[1]).length

                            if not (
                               points_original[0] == points_original[1] or
                               points_target[0] == points_target[1]
                               ):  # If any vector's length is zero

                                angle = vec_A.angle(vec_B) / pi
                            else:
                                angle = 0

                            # Set a range of acceptable variation in the connected edges
                            if dist_B > dist_A * 1.7 * self.join_stretch_factor or \
                               dist_B < dist_A / 2 / self.join_stretch_factor or \
                               angle >= 0.15 * self.join_stretch_factor:

                                merge_actual_vert = False
                                break
                    else:
                        merge_actual_vert = False
                except:
                    self.report({'WARNING'},
                        "Crosshatch set incorrectly")

                if merge_actual_vert:
                    coords = final_ob_duplicate.data.vertices[i].co
                    # To avoid problems when taking "-0.000" as a different value as "0.00"
                    for c in range(len(coords)):
                        if "%.3f" % coords[c] == "-0.00":
                            coords[c] = 0

                    comparison_coords = ["%.3f" % coords[0], "%.3f" % coords[1], "%.3f" % coords[2]]

                    if comparison_coords in main_object_verts_coords:
                        # Get the index of the vert with those coords in the main object
                        main_object_related_vert_idx = main_object_verts_coords.index(comparison_coords)

                        if self.main_object.data.vertices[main_object_related_vert_idx].select is True or \
                           self.main_object_selected_verts_count == 0:

                            ob_surface.data.vertices[i].co = final_ob_duplicate.data.vertices[i].co
                            ob_surface.data.vertices[i].select = True
                            crosshatch_verts_to_merge.append(i)

                            # Make sure the vert in the main object is selected,
                            # in case it wasn't selected and the "join crosshatch" option is active
                            self.main_object.data.vertices[main_object_related_vert_idx].select = True

        # Delete duplicated object
        bpy.ops.object.delete({"selected_objects": [final_ob_duplicate]})

        # Join crosshatched surface and main object
        bpy.ops.object.select_all('INVOKE_REGION_WIN', action='DESELECT')
        ob_surface.select_set(True)
        self.main_object.select_set(True)
        bpy.context.view_layer.objects.active = self.main_object

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

        bpy.ops.object.editmode_toggle('INVOKE_REGION_WIN')
        # Perform Remove doubles to merge verts
        if not (self.automatic_join is False and self.main_object_selected_verts_count == 0):
            bpy.ops.mesh.remove_doubles(threshold=0.0001)

        bpy.ops.mesh.select_all(action='DESELECT')

        # If the main object has modifiers, turn their "viewport view status"
        # to what it was before the forced deactivation above
        if len(self.main_object.modifiers) > 0:
            for m_idx in range(len(self.main_object.modifiers)):
                self.main_object.modifiers[m_idx].show_viewport = self.modifiers_prev_viewport_state[m_idx]

        self.update()

        return {'FINISHED'}

    def rectangular_surface(self, context):
        # Selected edges
        all_selected_edges_idx = []
        all_selected_verts = []
        all_verts_idx = []
        for ed in self.main_object.data.edges:
            if ed.select:
                all_selected_edges_idx.append(ed.index)

                # Selected vertices
                if not ed.vertices[0] in all_selected_verts:
                    all_selected_verts.append(self.main_object.data.vertices[ed.vertices[0]])
                if not ed.vertices[1] in all_selected_verts:
                    all_selected_verts.append(self.main_object.data.vertices[ed.vertices[1]])

                # All verts (both from each edge) to determine later
                # which are at the tips (those not repeated twice)
                all_verts_idx.append(ed.vertices[0])
                all_verts_idx.append(ed.vertices[1])

        # Identify the tips and "middle-vertex" that separates U from V, if there is one
        all_chains_tips_idx = []
        for v_idx in all_verts_idx:
            if all_verts_idx.count(v_idx) < 2:
                all_chains_tips_idx.append(v_idx)

        edges_connected_to_tips = []
        for ed in self.main_object.data.edges:
            if (ed.vertices[0] in all_chains_tips_idx or ed.vertices[1] in all_chains_tips_idx) and \
               not (ed.vertices[0] in all_verts_idx and ed.vertices[1] in all_verts_idx):

                edges_connected_to_tips.append(ed)

        # Check closed selections
        # List with groups of three verts, where the first element of the pair is
        # the unselected vert of a closed selection and the other two elements are the
        # selected neighbor verts (it will be useful to determine which selection chain
        # the unselected vert belongs to, and determine the "middle-vertex")
        single_unselected_verts_and_neighbors = []

        # To identify a "closed" selection (a selection that is a closed chain except
        # for one vertex) find the vertex in common that have the edges connected to tips.
        # If there is a vertex in common, that one is the unselected vert that closes
        # the selection or is a "middle-vertex"
        single_unselected_verts = []
        for ed in edges_connected_to_tips:
            for ed_b in edges_connected_to_tips:
                if ed != ed_b:
                    if ed.vertices[0] == ed_b.vertices[0] and \
                       not self.main_object.data.vertices[ed.vertices[0]].select and \
                       ed.vertices[0] not in single_unselected_verts:

                        # The second element is one of the tips of the selected
                        # vertices of the closed selection
                        single_unselected_verts_and_neighbors.append(
                                        [ed.vertices[0], ed.vertices[1], ed_b.vertices[1]]
                                        )
                        single_unselected_verts.append(ed.vertices[0])
                        break
                    elif ed.vertices[0] == ed_b.vertices[1] and \
                       not self.main_object.data.vertices[ed.vertices[0]].select and \
                       ed.vertices[0] not in single_unselected_verts:

                        single_unselected_verts_and_neighbors.append(
                                        [ed.vertices[0], ed.vertices[1], ed_b.vertices[0]]
                                        )
                        single_unselected_verts.append(ed.vertices[0])
                        break
                    elif ed.vertices[1] == ed_b.vertices[0] and \
                       not self.main_object.data.vertices[ed.vertices[1]].select and  \
                       ed.vertices[1] not in single_unselected_verts:

                        single_unselected_verts_and_neighbors.append(
                                              [ed.vertices[1], ed.vertices[0], ed_b.vertices[1]]
                                              )
                        single_unselected_verts.append(ed.vertices[1])
                        break
                    elif ed.vertices[1] == ed_b.vertices[1] and \
                       not self.main_object.data.vertices[ed.vertices[1]].select and \
                       ed.vertices[1] not in single_unselected_verts:

                        single_unselected_verts_and_neighbors.append(
                                              [ed.vertices[1], ed.vertices[0], ed_b.vertices[0]]
                                              )
                        single_unselected_verts.append(ed.vertices[1])
                        break

        middle_vertex_idx = None
        tips_to_discard_idx = []

        # Check if there is a "middle-vertex", and get its index
        for i in range(0, len(single_unselected_verts_and_neighbors)):
            actual_chain_verts = self.get_ordered_verts(
                                        self.main_object, all_selected_edges_idx,
                                        all_verts_idx, single_unselected_verts_and_neighbors[i][1],
                                        None, None
                                        )

            if single_unselected_verts_and_neighbors[i][2] != \
               actual_chain_verts[len(actual_chain_verts) - 1].index:

                middle_vertex_idx = single_unselected_verts_and_neighbors[i][0]
                tips_to_discard_idx.append(single_unselected_verts_and_neighbors[i][1])
                tips_to_discard_idx.append(single_unselected_verts_and_neighbors[i][2])

        # List with pairs of verts that belong to the tips of each selection chain (row)
        verts_tips_same_chain_idx = []
        if len(all_chains_tips_idx) >= 2:
            checked_v = []
            for i in range(0, len(all_chains_tips_idx)):
                if all_chains_tips_idx[i] not in checked_v:
                    v_chain = self.get_ordered_verts(
                                    self.main_object, all_selected_edges_idx,
                                    all_verts_idx, all_chains_tips_idx[i],
                                    middle_vertex_idx, None
                                    )

                    verts_tips_same_chain_idx.append([v_chain[0].index, v_chain[len(v_chain) - 1].index])

                    checked_v.append(v_chain[0].index)
                    checked_v.append(v_chain[len(v_chain) - 1].index)

        # Selection tips (vertices).
        verts_tips_parsed_idx = []
        if len(all_chains_tips_idx) >= 2:
            for spec_v_idx in all_chains_tips_idx:
                if (spec_v_idx not in tips_to_discard_idx):
                    verts_tips_parsed_idx.append(spec_v_idx)

        # Identify the type of selection made by the user
        if middle_vertex_idx is not None:
            # If there are 4 tips (two selection chains), and
            # there is only one single unselected vert (the middle vert)
            if len(all_chains_tips_idx) == 4 and len(single_unselected_verts_and_neighbors) == 1:
                selection_type = "TWO_CONNECTED"
            else:
                # The type of the selection was not identified, the script stops.
                self.report({'WARNING'}, "The selection isn't valid.")

                self.stopping_errors = True

                return{'CANCELLED'}
        else:
            if len(all_chains_tips_idx) == 2:    # If there are 2 tips
                selection_type = "SINGLE"
            elif len(all_chains_tips_idx) == 4:  # If there are 4 tips
                selection_type = "TWO_NOT_CONNECTED"
            elif len(all_chains_tips_idx) == 0:
                if len(self.main_splines.data.splines) > 1:
                    selection_type = "NO_SELECTION"
                else:
                    # If the selection was not identified and there is only one stroke,
                    # there's no possibility to build a surface, so the script is interrupted
                    self.report({'WARNING'}, "The selection isn't valid.")

                    self.stopping_errors = True

                    return{'CANCELLED'}
            else:
                # The type of the selection was not identified, the script stops
                self.report({'WARNING'}, "The selection isn't valid.")

                self.stopping_errors = True

                return{'CANCELLED'}

        # If the selection type is TWO_NOT_CONNECTED and there is only one stroke, stop the script
        if selection_type == "TWO_NOT_CONNECTED" and len(self.main_splines.data.splines) == 1:
            self.report({'WARNING'},
                        "At least two strokes are needed when there are two not connected selections")

            self.stopping_errors = True

            return{'CANCELLED'}

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

        bpy.ops.object.select_all('INVOKE_REGION_WIN', action='DESELECT')
        self.main_splines.select_set(True)
        bpy.context.view_layer.objects.active = self.main_splines

        # Enter editmode for the new curve (converted from grease pencil strokes), to smooth it out
        bpy.ops.object.editmode_toggle('INVOKE_REGION_WIN')
        bpy.ops.curve.smooth('INVOKE_REGION_WIN')
        bpy.ops.curve.smooth('INVOKE_REGION_WIN')
        bpy.ops.curve.smooth('INVOKE_REGION_WIN')
        bpy.ops.curve.smooth('INVOKE_REGION_WIN')
        bpy.ops.curve.smooth('INVOKE_REGION_WIN')
        bpy.ops.curve.smooth('INVOKE_REGION_WIN')
        bpy.ops.object.editmode_toggle('INVOKE_REGION_WIN')

        self.selection_U_exists = False
        self.selection_U2_exists = False
        self.selection_V_exists = False
        self.selection_V2_exists = False

        self.selection_U_is_closed = False
        self.selection_U2_is_closed = False
        self.selection_V_is_closed = False
        self.selection_V2_is_closed = False

        # Define what vertices are at the tips of each selection and are not the middle-vertex
        if selection_type == "TWO_CONNECTED":
            self.selection_U_exists = True
            self.selection_V_exists = True

            closing_vert_U_idx = None
            closing_vert_V_idx = None
            closing_vert_U2_idx = None
            closing_vert_V2_idx = None

            # Determine which selection is Selection-U and which is Selection-V
            points_A = []
            points_B = []
            points_first_stroke_tips = []

            points_A.append(
                self.main_object.matrix_world @ self.main_object.data.vertices[verts_tips_parsed_idx[0]].co
                )
            points_A.append(
                self.main_object.matrix_world @ self.main_object.data.vertices[middle_vertex_idx].co
                )