This post is inspired by a happy coincidence, I was doing p2 of day 14 of advent of code for 2017 and I had never formerly encountered the flood fill algorithm. This confluence of events inspired me to come up with a new way to solve this classic problem.
You are given a boolean matrix and told that points are neighbors if they are either next to each other horizontally or vertically, but not diagonally. If point (a) and point (b) are neighbors and point (c) is a neighbor of point (b) then (a),(b) and (c) are in the same neighborhood. A point with no neighbors is a neighborhood of 1. You are asked how many such neighborhoods exist.
Shamelessly borrowing the picture from adventofcode:
11.2.3..-->
.1.2.3.4
....5.6.
7.8.55.9
.88.5...
88..5...
.8...8..
8888888.-->
| |
V V
You can see that this piece of the matrix, has 9 neighborhoods. Points that are off in the boolean matrix are simply ignored.
This problem is a classic because it is apparently how that bucket fill tool works in every paint program. Essentially you keep filling all the pixels that are neighbors until you are done.
Now the classical way to solve this problem is flood fill. Essentially the problem is broken into two pieces:
- An outer loop that visits every point
- An inner loop that looks for neighbors using breadth first search or depth first search
Here I will annotate the Q code that achieves this result.
ff:{[m] /this is flood fill function
/ returns a table with all the distinct neighborhoods
/m is a matrix
[m]
/c is the number of rows
/columns this matrix has, we assume it's square
c:count m;
/we flatten m; m is now a list of length c*c
m:raze m;
/we create t, a table that is empty
/and has columns id, and n
/id will be the id of the neighborhood
/n will be a list of neighbors in that hood
t:([id:()]n:());
/gni is a helper function that get neighbor indexes
/gni will accept the index as an offset
/ from the flattened matrix
/and will return coordinates as flat offsets.
/internally it needs to represent
/the coordinates as a pair,
/ we use standard div and mod to represent the pair
gni: { [x;y]
/x is the size of the matrix it will be fixed to c
/y is the point
/s whose neighbors we want to find
/ny are the new col coordinates
/ found by dividing the point by c
/ we add the offsets
/ but makes sure we are within boundaries of matrix
/ min 0 and max c-1
ny:0|(x-1)&0 -1 1 0+/:y div x;
/repeat this process to calculate row coordinates
nx:0|(x-1)&-1 0 0 1+/:y mod x;
/get back original flattened index
/by multiplying ny by c and adding nx
/flatten list and only return distinct points
distinct raze nx+x*ny
/project function onto c, since it will be constant
/ for the duration of the parent function
}[c];
/f is our breadth first search function
/we start by adding i, which is our first point,
/we then find all the gni of i,
/we remove any gni[i] that are already present in i
/and intersect the new points that are on (=1)
/by intersecting with j
/we project this function on gni,
/ which is called gni inside as well.
/we also project on the points that are on
/(are equal to 1) aka (where m)
/ and which are named j inside the function)
f:{[gni;i;j](),i,j inter gni[i] except i}[gni;;where m];
/repeating f over and over will yield
/ all the neighbors of a point that are on (=1)
/now we get to the outer loop,
/which will check
/if we have seen the current point (p)
/if we have seen it, we ignore it,
/and return our current table (t)
/otherwise, we call f on that point and return
/the point and all of it's neighbors
/this will end up with our new table t,
/we only call this function on the points
/that are turned on,
/which means we can skip checking if the point is on
/we will repeat this for all the points equal to 1
t:{[f;m;t;p]$[not p in exec n from t;
t upsert (p;f over p);
t]}[f;n]/[t;where m];
}
The entire code without comments looks like this, and is pretty short.
ff:{[m]
c:count m;m:raze m; t:([id:()]n:());
gni:{distinct raze (0|(x-1)&-1 0 0 1+/:y mod x)+x* 0|(x-1)&0 -1 1 0+/:y div x}[c];
f:{[gni;i;j](),i,k where j k:gni[i] except i}[gni;;m];
t:{[f;m;t;p]$[not p in raze exec n from t;t upsert (p;f over p);t]}[f;n]/[t;where m]
}With that out of the way, I now want to explain an alternate way of achieving something essentially identical.
A common technique in math and engineering fields is if you can’t solve a problem, find a similar problem that is easier that you can solve and solve that first.
So I pretended that all I needed to do was to find runs, which is a set of consecutive 1s in a row. I figured if I do this for all the rows and all the columns, I’ll have a starter set of neighborhoods and I can figure how to merge them later.
Finding a run is really simple:
- set the neighborhood id to 0
- start at the beginning of first row
- look at each element,
- if you see a 1 add element to the list that corresponds to your current neighborhood id
- if you see a zero increment the neighborhood id
- if next row, go to the next row, increment neighborhood id, otherwise go to step 8
- go to step 3
- remove neighborhood ids with no elements in them.
In q, you can take a bit of a shortcut by asking each row, for the indexes with a 1 in them. So for example:
q)where 1 0 1 1 0
0 2 3
q) deltas where 1 0 1 1 0
0 2 1
One thing we notice is that consecutive differences of 1 represent runs and those with any other number represent cut points, and if we cut at those points, we would get (,0; 2 3) which is a list that contains two lists, each inner list is a run in that row. We can repeat this for the columns and create all the runs that exist row wise or column wise.
Now we just need some logic to merge these runs. Here is an annotated version of the code that implements this method of finding all the neighborhoods.
rr:{[m] /This will return all the points in the matrix that are on and their assigned neighborhood id
/m is the square matrix
c:count m; /c is the length of a side of the matrix
m:raze m; /m is now the flattened list that corresponds to the matrix
/runs is a function that will return a list of lists,
/each inner list corresponds to points that are adjacent
/Though the mechanism to calculate adjacency is in terms of y
/It will return the Xs that correspond to those Ys.
/eg if we have a row 0 1 1 0 1 0 1 1, and the points are labeled `a`b`c`d`e`f`g`h
/runs[0 1 1 0 1 0 1 1;`a`b`c`d`e`f`g`h] => (`b`c;,`e;`g`h)
runs:{(where 1<>-2-': y) _ x};
/now we want to apply runs to every row and column
/n will be our initial list of neighborhoods,
/ each neighborhood is a list that will contain adjacent elements either by row or by column
/ elements are given a unique id
/i will be our unique id , it is the index of the point that is turned on (=1)
i:where m;
/j groups all 'i's that are in the same row, by using integer division
/so j is a dictionary whose keys are the row index and whose values are the 'i's
/we need to do this grouping to ensure that only points in the same row are considered adjacent
/ otherwise we might believe that a diagonal matrix has runs (because we only see the column indexes)
j:i group i div c;
/since runs returns x, but uses 'y's to do the splitting, we provide j mod c as y
/j mod c corresponds to the column location of each i
/we then apply runs on (')each j and j mod c.
/finally we flatten the result with a raze so that we just a have the list of row neighborhoods
n:raze runs'[j;j mod c];
/now we will repeat with the columns
/we want to permute the list i so that the elements appear in column order
/we want to preserve the labels we have in i
/so we generate a list that is the length of m (c*c),
/we reshape (#) to be a square matrix, and transpose it along the diagonal (flip)
/we flatten(raze) it back out and intersect it with i,
/this preserves the column order we achieved from doing the flip,
/but only keeps the elements that are in i
/k are the elements of i permuted along columns
/eg if we have a matrix
// 0 0 1
// 1 0 0
// 0 1 0
// i would be 2 3 7
// and k would be 3 7 2
k:(raze flip (c;c)#til c*c) inter i);
/now we will create the j again which is our grouping
/but this time according to the column index so we will group by the integer remainder (mod)
/this will separate our indexes into columns
j:k group (k mod c);
/then we will apply our runs function on each column
/ we use the div to get the row index of each element,
/ so that runs can tell if two points are vertically adjacent
/ we flatten the result and append it to n
n,:raze runs'[j;j div c];
/ at this point we have a list of all possible neighborhoods, but there are overlaps
/ so we create a table t, which will assign an id for each neighborhood contender
/ ungroup flattens the table so that every element belongs to a neighborhood id
// here is a sample of what the table looks like before and after ungroup
/ before ungroup:
/ n id
/ ------------
/ 0 1 0
/ ,3 1
/ ,5 2
/ 8 9 10 11 3
/after ungroup
/ id n
/ -----
/ 0 0
/ 0 1
/ 1 3
/ 2 5
/ 3 8
/ 3 9
/ 3 10
/ 3 11
t:ungroup `id xkey update id:i from ([]n);
/ f is our merging function.
/ It will merge neighborhoods together,
/ by giving them the same id if there are points in common.
/ First we will update/create a column p which is the minimum id for a given element (n)
/ Then we will update all of the ids the with minimum p and rename that the new id
/ If two neighborhoods share a point in common, that point has two neighborhood id's
/ we will assign it primary id (p) which is the min of those two ids
/ we then group by id and assign the id to be the min of the primary ids (p)
f:{[t]update id:min p by id from update p:min id by n from t};
/On each iteration of this function we will merge some neighborhoods,
/ we repeat until we can't merge anymore
f over t}
The uncommented code is much shorter and looks like this:
rr2:{[m]
c:count m;m:raze m; runs:{(where 1<>-2-': y) _ x};
n:raze runs'[j;(j:i group (i:where m) div c) mod c];
n,:raze runs'[j;(j:k group (k:(raze flip (c;c)#til c*c) inter i) mod c) div c];
t:ungroup `id xkey update id:i from ([]n);
f:{[t]update id:min p by id from update p:min id by n from t};
f over t}
At this point we are done and we have solved the puzzle in two slightly different ways. Looking at the performance of both methods, reveals that they are both O(n) essentially linear in the number of elements. However, the row method is actually faster, because everything is done as nice straight rows, so there is not much jumping around, in other words it is more cpu cache friendly. On my late 2011 Mac: for a 128 by 128 matrix
\t rr m
176 milliseconds
\t ff m
261 milliseconds
So it is about 66% faster
As a bonus, my father suggested that I expand this code into n dimensions. This is rather easy to do and I include the code here without too many comments. The trick is taking advantage of KDB’s vs (vector from scalar) function that can rewrite an integer as a number in another base:
so for example 3 vs 8 => (2;2) which are the coordinates of the point 8 in a 3 by 3 matrix. The point 25 is 2 2 1, in other words it’s in the last layer of a 3 dimensional cube of size 3, last column and first row.
This means we can easily go back and forth between integer representation of the matrix and coordinate representation. Everything else is pretty much the same.
ffN:{[m]
c:count m;r:1; while[(count m:raze m)<>count[m];r+:1]; t:([id:()]n:());
gniN:{[r;c;n]a:neg[a],a:a where 1=sum each a:(cross/)(r;2)#01b;
distinct raze -1 _ flip c sv flip (c-1)&0|a +\: c vs n,-1+prd r#c}[r;c];
f:{[gni;i;j](),i,j inter gni[i] except i}[gniN;;where m];
ps:where m;
while[count ps;t:t upsert (first ps;k:f over first ps); ps:ps except k];
/functional version is slower
/t:{[f;t;j]q+::1;$[not j in raze exec n from t;t upsert (j;f over j);t]}[f]/[t;where m]
t}
rrN:{[m]
c:count m;r:1; while[(count m:raze m)<>count[m];r+:1]; runs:{(where 1<>-2-': y) _ x};
j:(c sv (til r) rotate\: c vs til prd r#c) inter\: where m;
g:til[r] except/: d:reverse til r;
n:raze raze each runs''[j@'k;@'[v;d]@'k:('[group;flip]) each @'[v:flip c vs j;g]];
t:ungroup `id xkey update id:i from ([]n);
f:{[t]update id:min p by id from update p:min id by n from t};
f over t}
iabdb 