How can we find the largest contiguous region of a graph with two different "ID"s?

Question

I've recently learned about the Flood-Fill Algorithm, an algorithm that can take a graph and assign each node a component number in O(N) time.

For example, a common problem that can be solved efficiently with the Flood-Fill Algorithm would be to find the largest region in a N*N board, where every node in the region is adjacent to another node with the same ID either directly up, down, to the left, or to the right.

In this board, the largest regions would both be of size 3, made up of all 1s and all 9s respectively.

However, I recently started wondering if we could extend this problem; specifically, if we could find the largest region in a graph such that every node in the region is adjacent to another node with two possible IDs. In the above board, the largest such region is made up of 1s and 9s, and has a size of 7.

Here was my thought process in trying to solve this problem:

Thought 1: O(N^4) Algorithm

We can solve this in O(N^4) time using a basic flood-fill algorithm. We do this by testing all O(N^2) pairs of horizontally or vertically adjacent squares. For every pair of squares, if they have different IDs, then we run a flood-fill from one of the two squares.

Then, by modifying the flood-fill algorithm so that it travels to squares with one of the two possible IDs, we can test each pair in O(N^2) time --> O(N^2) pairs * O(N^2) flood fill per pair = O(N^4) algorithm.

Then, I had an insight: An Possibly O(N^2) Algorithm

First, we run a regular flood-fill through the board and separate the board into a "component graph" (where each component in the original graph is reduced to a single node).

Now, we do a flood-fill through the edges of the component graph instead of the nodes. We mark each edge with a pair of integers signifying the two IDs inside the two components which it connects, before flood-filling through the edges as if they themselves were nodes.

I believe that this, if implemented correctly, would result in a O(N^2) algorithm, because an upper bound for the number of edges in a N*N board is 4*N*N.

Now, my question is, is my thought process logically sound? If not, can somebody suggest another algorithm to solve this problem?

Answer 1

Here is the algorithm that I wrote to solve your problem. It expands on your idea to flood-fill through the edges (great idea, by the way) and is able to output the correct answer for a 250*250 grid in less than 300ms, with less than 30 megabytes of memory allocated.

Here is the problem that I managed to find online that matches your question exactly, and it is also where I tested the validity of my algorithm: USACO Problem

Note that the USACO Problem requires us to find the largest single-id component before finding the largest double-id component. In my algorithm, the first step is actually necessary in order to reduce the whole board into a component graph.

Here's my commented C++ Code:

#include <iostream>
#include <fstream>
#include <cmath>
#include <algorithm>
#include <vector>
#include <unordered_set>

using namespace std;

// board to hold square ids and comp[][] to mark component numbers
vector <vector<int>> board, comp;
vector <int> comp_size = {-1}; // size of those components
vector <int> comp_id = {-1}; // id contained within those components
vector <unordered_set <int>> adj = {{}}; // component graph adjacency list
vector <bool> visited; // component graph visited array

void dfs(int x, int y, int N, int id, int curr_comp){
    if(x < 0 || x >= N || y < 0 || y >= N){return;}
    else if(board[x][y] != id){
        if(comp[x][y] == 0){return;}

        // construct component graph adjacency list during the first flood-fill
        adj[comp[x][y]].insert(curr_comp);
        adj[curr_comp].insert(comp[x][y]);
        // this is why we use an unordered_set: it automatically eliminates
        // collisions

        return;
    }
    else if(comp[x][y]){return;}

    ++comp_size[curr_comp];
    comp[x][y] = curr_comp;

    dfs(x-1, y, N, id, curr_comp);
    dfs(x+1, y, N, id, curr_comp);
    dfs(x, y-1, N, id, curr_comp);
    dfs(x, y+1, N, id, curr_comp);
}

void dfs2(int curr, int id1, int id2, int &size){
    visited[curr] = true;

    // recurse from all valid and adjacent components to curr
    vector <int> to_erase;
    for(int item : adj[curr]){
        if(visited[item]){continue;}
        if(comp_id[item] == id1 || comp_id[item] == id2){
            to_erase.push_back(item);
            size += comp_size[item];
            dfs2(item, id1, id2, size);
        }
    }

    // we erase all edges connecting the current component AT THE SAME TIME to
    // prevent std::unordered_set iterators from being invalidated, which would
    // happen if we erased items as we iterated through adj[curr]
    for(int item : to_erase){
        adj[curr].erase(item);
        adj[item].erase(curr);
    }

    return;
}

int main()
{
    ifstream fin("multimoo.in");
    ofstream fout("multimoo.out");

    int N;
    fin >> N;

    board = vector <vector<int>> (N, vector <int> (N));
    for(int i = 0; i < N; ++i){
        for(int j = 0; j < N; ++j){
            fin >> board[i][j];
        }
    }
    // Input Done

    comp = vector <vector<int>> (N, vector <int> (N, 0)); // note that comp[i][j] = 0 means not visited yet

    // regular flood-fill through all the nodes
    for(int i = 0, curr_comp = 1; i < N; ++i){
        for(int j = 0; j < N; ++j){
            if(comp[i][j]){continue;}

            // add information about the current component
            comp_size.push_back(0);
            comp_id.push_back(board[i][j]);
            adj.push_back({});

            dfs(i, j, N, board[i][j], curr_comp++);
        }
    }

    fout << *max_element(comp_size.begin(), comp_size.end()) << endl;

    int ANS = 0;
    for(unsigned int i = 1; i < comp_size.size(); ++i){
        // no range-for loop here as we erase elements while iterating, which
        // may invalidate unordered_set iterators; instead, we use a while-loop
        while(!adj[i].empty()){
            int size = comp_size[i], curr = *(adj[i].begin());
            visited = vector <bool> (comp_size.size(), false); // reset visited
            dfs2(i, comp_id[i], comp_id[curr], size);
            ANS = max(ANS, size);
        }
    }

    fout << ANS << endl;

    return 0;
}

As for the time complexity, I personally am not very sure. If somebody could help analyze this algorithm to determine its complexity, I'd greatly appreciate it!

Answer 2

Your algorithm works...

As far as I can tell, flood filling over your induced graph indeed gives all possible components, after which it's simple to find the largest one.

...but I'm not sure about the runtime

You correctly say that there are O(N^2) edges in the original graph, and therefore O(N^2) nodes in the induced graph. However, these nodes are no longer guaranteed to be in a nice grid, which may leave more than O(N^2) induced edges.

For example, consider the large "1-block" in your example. This block has 6 edges, which will give a complete graph with 6 vertices, as all these edges-turned-vertices are connected. This may give you an induced graph with more than O(N^2) edges, making it impossible to find components in O(N^2) time.

Therefore, I believe that the algorithm will not run in O(N^2), but I'm unsure of the actual runtime, as it will depend on what exactly the algorithm does at this point. The question only notes flood fill, but I think it had not imagined this situation.

Consider the following 9x9 grid:

232323232
311111113
212313212
313212313
212313212
313212313
212313212
313212313
212313212

The idea is simple: it's a single large component designed to border as many small components as possible. The induced graph here would be a single almost-complete graph with O(N^2) vertices and O(N^4) edges. Alternatively, if we only link the (1,2) edges with other (1,2) edges, and similar for (1,3) edges and other (1,3) edges, we will have a slightly less-connected graph, but it would still consist of two components with each O(N^4) edges, albeit with a lower constant.

Therefore, creating this graph would take at least O(N^4) time, as would traversing it. This is the time I would argue that the algorithm takes, but I cannot prove that there are no possible optimizations that improve upon this.

Answer 3

We could achieve the optimal O(N^2) complexity by smartly doing our DFS from each pivot component.

Explanation:

1. First create the set of same-valued components and relationship-map to their neighbours
2. Notice that for the flood-zone we are looking for atmost 2 distinct values.
3. Lets say we consider the point (i, j) and look at its neighbours
4. For each 2-value-pair, say [v_ij, v_neighbour] => do a BFS from this (i,j) pivot point while only collecting nodes such that node-value is one of [v_ij, v_neighbour]
5. Notice that each component is visited only constant times for BFS (we ensure that by deleting the reverse-edge from child->parent while doing BFS).
6. Because of (5), our complexity remains O(N^2)

Working code in Python:

from queue import Queue


class Comp:
    def __init__(self, point, value):
        self.members = {point}
        self.value = value
        self.neighbours = set()
        self.pivot = point

    def can_add_member(self, value):
        return value == self.value

    def add_member(self, point):
        self.members.add(point)

    def add_neighbour(self, neighbour_comp):
        self.neighbours.add(neighbour_comp)

    def __str__(self):
        return f'[M:%d, V:%d, N:%d]' % (len(self.members), self.value, len(self.neighbours))


def find_largest_flood_region(D):
    point_to_comp_map = {}
    N, M = len(D), len(D[0])

    # Step-1: Create same-value connected-components:
    for x in range(N):
        for y in range(M):
            if (x, y) in point_to_comp_map:
                continue
            comp_xy = Comp((x, y), D[x][y])
            point_to_comp_map[(x, y)] = comp_xy
            pq = Queue()
            pq.put((x, y))
            while pq.qsize() > 0:
                i, j = pq.get()
                for l, m in [(i-1, j), (i+1, j), (i, j-1), (i, j+1)]:
                    if 0 <= l < N and 0 <= m < M and (l, m) not in point_to_comp_map and D[l][m] == D[x][y]:
                        comp_xy.add_member((l, m))
                        point_to_comp_map[(l, m)] = comp_xy
                        pq.put((l, m))

    # Step-2: Create the relationship-map between the components created above
    for x in range(N):
        for y in range(M):
            comp_xy: Comp = point_to_comp_map[(x, y)]
            for i, j in [(x-1, y), (x+1, y), (x, y-1), (x, y+1)]:
                if 0 <= i < N and 0 <= j < M and D[i][j] != D[x][y]:
                    comp_ij: Comp = point_to_comp_map[(i, j)]
                    comp_xy.add_neighbour(comp_ij)
                    comp_ij.add_neighbour(comp_xy)

    # Do BFS one by one on each unique component:
    unique_comps = set(point_to_comp_map.values())
    max_region = 0
    for comp in unique_comps:
        potential_values = set([neigh_comp.value for neigh_comp in comp.neighbours])
        for value in potential_values:
            value_set = {value, comp.value}
            region_value = 0
            pq = Queue()
            pq.put(comp)
            while pq.qsize() > 0:
                comp_xy: Comp = pq.get()
                region_value += len(comp_xy.members)
                for ncomp in comp_xy.neighbours:
                    if ncomp.value in value_set:
                        if comp_xy in ncomp.neighbours:
                            ncomp.neighbours.remove(comp_xy)
                        pq.put(ncomp)
            max_region = max(max_region, region_value)

    return max_region


D = [
    [9,2,7,9],
    [1,1,9,9],
    [3,1,4,5],
    [3,5,6,6]
]

print(find_largest_flood_region(D))

Output:

7

Answer 4

We can show that solving this in O(n), where n is the number of elements in the matrix, is possible with two passes of a flood-fill union-find routine without a depth-first search.

Given

9 2 7 9
1 1 9 9
3 1 4 5
3 5 6 6

after we label with flood fill, we have:

A B C D
E E D D
F E G H
F I J J

Now that we know each component's size, we can restrict each cell to testing its best connection to a different component left or up. We only need to check one field in a map on the component, the one pointing to the same number, potentially create a new component of reference, or merge two.

In the following example, we'll label components with more than one value with two letters unrelated to their original components. Each cell visited can generate at most two new components and updates at most two components so the complexity remains O(n).

Iterating left to right, top to bottom:

A0: {⊥: 1, 2: AA, 1: DD}

B0: {⊥: 1, 9: AA, 7: BB, 1: EE}
  AA = {size: 2}

C0: {⊥: 1, 2: BB, 9: CC}
  BB = {size: 2}

D0: {⊥: 3, 7: CC}
  CC = {size: 4}

E0: {⊥: 3, 9: DD}
  DD = {size: 4}

E1: {⊥: 3, 9: DD, 2: EE}
  EE = {size: 4}

D1: {⊥: 3, 7: CC, 1: DD}
  DD updates to size 7

D2: {⊥: 3, 7: CC, 1: DD}

F0: {⊥: 2, 1: FF}
  FF = {size: 5}
  
... etc.