def cover_edges(cs):
    pts = sorted(cs, key=lambda p: p[0], reverse=True)   # decreasing x
    # maintain processed points sorted by y as list of (y,x) tuples
    import bisect
    S = []            # sorted by (y, x)
    edges = []
    for x, y in pts:
        # find first processed point with y' > y
        idx = bisect.bisect_right(S, (y, float('inf')))
        min_x = float('inf')
        # scan upward in increasing y
        i = idx
        while i < len(S):
            yq, xq = S[i]
            if xq < min_x:
                edges.append(((x, y), (xq, yq)))
                min_x = xq
            i += 1
        # insert current point into S maintaining sort-by-y
        bisect.insort(S, (y, x))
    return edges

cs = [[1,3],[2,2],[3,1],[4,4]]
print(cover_edges(cs))