Im trying to re-implement the 3d seam carving from this paper in python: https://people.csail.mit.edu/mrub/papers/vidret.pdf

Specifically I have a 3d array called mask_data that masks out areas of the 3d volume I want to be prioritized, aka not removed. For testing I just use an array of 0’s with a 20% chance any of the places become a 1. The paper shows how to theoretically construct a graph that will allow the min cut to be a decent seam to remove, and the following function does that. The issue is that even with a 128^3 graph it takes 13 minutes to construct the graph on my M2 max macbook. I tried parallelising the operation but I only succeeded in making it take longer.

Any ideas/implementations that would allow for faster creation of the graph structure? I am using the graph_tool library that implements its functions in C++ and calculates the max flow min cut quickly in the next step so id like to continue using it.

import graph_tool.all as gt
def create_directed_energy_graph_from_mask(mask_data, direction='left', large_weight=1e8):
    z,y,x = mask_data.shape  # Dimensions of the 3D mask array
    g = gt.Graph(directed=True)
    weight_prop = g.new_edge_property("int")  # Edge property for weights

    # Function to get linear index from 3D coordinates, assuming C-style row-major order
    index = lambda i, j, k: i * y * x + j * x + k

    # Add vertices
    num_vertices = z * y * x
    g.add_vertex(num_vertices)
    vertices = list(g.vertices())  # Get the list of vertices

    # Define neighbor offsets based on directionality
    # Source to sink propagation direction - positive axes directions
    directions = {
        'left': [(0, 0, 1)],  # propagate right
        'right': [(0, 0, -1)],  # propagate left
        'top': [(0, 1, 0)],  # propagate downwards
        'bottom': [(0, -1, 0)],  # propagate upwards
        'front': [(1, 0, 0)],  # propagate back
        'back': [(-1, 0, 0)]  # propagate front
    }

    neighbors = directions[direction]
    print(x,y,z)
    
    for i in range(z):
        for j in range(y):
            for k in range(x):
                current_index = index(i, j, k)
                # print(current_index, len(vertices))
                current_vertex = vertices[current_index]

                # Check each neighbor direction for valid connections
                for di, dj, dk in neighbors:
                    ni, nj, nk = i + di, j + dj, k + dk
                    if 0 <= ni < z and 0 <= nj < y and 0 <= nk < x:
                        neighbor_index = index(ni, nj, nk)
                        neighbor_vertex = vertices[neighbor_index]

                        # print(current_index, neighbor_index)

                        # Determine edge weight
                        weight = 10 if mask_data[ni, nj, nk] != 0 or mask_data[i, j, k] != 0 else 1

                        # Add edge and assign weight
                        e = g.add_edge(current_vertex, neighbor_vertex) #forward edge with energy value
                        e2 = g.add_edge(neighbor_vertex, current_vertex) #backward edge with large energy value
                        weight_prop[e] = weight
                        weight_prop[e2] = large_weight

                # Add each diagonal backwards neighbor inf edge, ie x-1, y-1 and x-1, y+1 for YX plane
                if k > 0 and j > 0:
                    neighbor_index = index(i, j-1, k-1)
                    neighbor_vertex = vertices[neighbor_index]
                    # print(current_index, neighbor_index)
                    # print(i,j,k, " y-1, x-1 ", i, j-1, k-1) 
                    e = g.add_edge(current_vertex, neighbor_vertex)
                    weight_prop[e] = large_weight+1
                if k > 0 and j < y-1:
                    neighbor_index = index(i, j+1, k-1)
                    neighbor_vertex = vertices[neighbor_index]
                    # print(current_index, neighbor_index)
                    # print(i,j,k, "  y+1 x-1 ", i, j+1, k-1)
                    e = g.add_edge(current_vertex, neighbor_vertex)
                    weight_prop[e] = large_weight+1

                # print(i,j,k, current_index)
                    
                # Add each diagonal backwards neighbor inf edge for ik plane
                if k > 0 and i > 0:
                    neighbor_index = index(i-1, j, k-1)
                    neighbor_vertex = vertices[neighbor_index]
                    # print(current_index, neighbor_index)
                    e = g.add_edge(current_vertex, neighbor_vertex)
                    weight_prop[e] = large_weight+2
                if k > 0 and i < z-1:
                    neighbor_index = index(i+1, j, k-1)
                    neighbor_vertex = vertices[neighbor_index]
                    # print(current_index, neighbor_index)
                    e = g.add_edge(current_vertex, neighbor_vertex)
                    weight_prop[e] = large_weight+2

    g.edge_properties["weight"] = weight_prop
    return g, weight_prop

The paper does mention multi-resolution banding as a speed up, but all the other parts of the pipeline run quickly, its just creating the graph that is slow. I may have to figure out what they mean by that though if I cant optimise this code. Also note that this version only works for the ‘left’ direction as of now.

I tried parallelising the code by splitting the mask data array and computing subgraphs, expected it to run faster, but it ran slower, perhaps due to the graph recombination operation being expensive.