# Bridge the Gaps

Given a black and white image with a white background and a set of black dots, paint a set of white pixels red, such that there is a path between each pair of black pixels.

### Details

• A path is a set of connected pixels (8-neighbourhood connectivity). Black pixels can be used as part of the paths. The goal is trying to minimize the set of red pixels under the above conditions, and outputting a corresponding image.

• You don't have to find the optimal solution.

• A trivial and at the same time worst solution is just painting all the white pixels red.

• Example (Pixels are enlarged for visibility): ## Details

• Given a pixel image (in any suitable format) return another image with the dots connected as specified above, as well as a integer indicating how many red pixel were used.
• The Score is the product of (1+the number of red pixels) for each of the 14 testcases.
• The goal is having the lowest score.

## Testcases

The 14 testcases are shown below. A python program to verify the connectedness of the outputs can be found here.              ## Meta

Thanks to @Veskah, @Fatalize, @wizzwizz4 and @trichoplax for the various suggestions.

• Good challenge; I like ones with different and creative scoring schemes. I assume the program needs to work on an arbitrary image, not just these 14 specific examples? If so, can we assume a reasonable max size, like 512x512 per the Mona Lisa image, or 1024x1024? Aug 5 '19 at 14:59
• Thanks for the feedback! Yes you can assume a maximal size (an also minimal size if necessary), as long as all the 14 examples can be processed. Aug 5 '19 at 15:01
• how do i convert png to ascii or json or something else easy to parse?
– ngn
Aug 5 '19 at 18:26
• Do you have to be able to compute your own score? A program which tries every possible combination of white pixels to paint red, and sees which subset has the fewest red pixels while connecting all the black pixels would have the best possible score, but it would be so slow that it would take longer than the lifetime of the universe to actually compute that score. Aug 6 '19 at 0:01
• @ngn Open in GIMP, save as netpbm format. Aug 7 '19 at 12:52

# Python, 2.62 * 10^40

This algorithm just floodfills (BFS) the plane starting from the black parts of the image, where for each new pixel we record what black part it was flooded from. As soon as we have two neighbouring pixels with different black parts as ancestors, we basically merge these two black parts by joining them through the ancestors of the two neighbours we just found. In theory this could be implemented in O(#pixels), but to keep the amount of code at an acceptable level this implementation is slightly worse. ### Output

import numpy as np
from scipy import ndimage
import imageio
from collections import deque

for k in range(1, 15):
fname=str(k).zfill(2) +'.png'
print("processing ", fname)

print(img.shape)

# determine non_white part
white = np.logical_and(np.logical_and(img[:,:,0] == 255, img[:,:,1] == 255), img[:,:,2] == 255)
non_white = np.logical_not(white)

# find connected components of non-white part
neighbourhood = np.ones((3,3))
labeled, nr_objects = ndimage.label(non_white, neighbourhood)

# print result
print("number of separate objects is {}".format(nr_objects))

# start flood filling algorithm
ind = np.nonzero(labeled)
front = deque(zip(ind,ind))

membership = np.copy(labeled)
is_merge_point = np.zeros_like(labeled) > 0
parent = np.zeros((2,) + labeled.shape) #find ancestor of each pixel
is_seed = labeled > 0
size_i, size_j = labeled.shape
# flood from every seed
while front: #while we have unexplored pixels
point = front.popleft()
# check neighbours:
for (di,dj) in [(-1,-1),(-1,0),(-1,1),(0,-1),(0,1),(1,-1),(1,0),(1,1)]:
current = membership[point, point]
new_i, new_j = point+di, point+dj
if 0 <= new_i < size_i and 0 <= new_j < size_j:
value = membership[new_i, new_j]
if value == 0:
membership[new_i, new_j] = current
front.append((new_i, new_j))
parent[:, new_i, new_j] = point
elif value != current: #MERGE!
is_merge_point[point, point] = True
is_merge_point[new_i, new_j] = True
membership[np.logical_or(membership == value, membership == current)] = min(value, current)

# trace back from every merger
ind = np.nonzero(is_merge_point)
merge_points = deque(zip(ind.astype(np.int),ind.astype(np.int)))
for point in merge_points:
next_p = point
while not is_seed[next_p, next_p]:
is_merge_point[next_p, next_p] = True
next_p = parent[:, next_p, next_p].astype(np.int)

img_backup = np.copy(img)
img[:,:,0][is_merge_point] = 255 * img_backup[:,:,0]
img[:,:,1][is_merge_point] = 0   * img_backup[:,:,1]
img[:,:,2][is_merge_point] = 0   * img_backup[:,:,2]

#compute number of new points
n_red_points = (img[:,:,0] != img[:,:,1]).sum()
print("#red points:", n_red_points)

# plot: each component should have separate color
imageio.imwrite("./out_images/"+fname, np.array(img))



### Score

(1+183)*(1+142)*(1+244)*(1+42)*(1+1382)*(1+2)*(1+104)*(1+7936)*(1+26)*(1+38562)*(1+42956)*(1+6939)*(1+8882)*(1+9916)
= 26208700066468930789809050445560539404000
= 2.62 * 10^40

• --This, I believe, is optimal. Well done.-- Okay, this isn't optimal. I don't understand why not. Aug 7 '19 at 19:31
• @wizzwizz4 Look at the easy case of the four corners of a square: The optimal solution would be an X. While in theory my algorithm could find this solution, it is very unlikely. It is a lot more likely that it finds a solution with three paths each connecting two points. Aug 7 '19 at 20:37
• @wizzwizz4 Yep, zoom in on the wikipedia text example, and you'll see tons of little places where a different connecting path would have saved a red pixel or two; they'll add up. Aug 7 '19 at 21:10
• But this seems like soap-bubbles on pegs, which is a legitimate solution to the Steiner tree problem. Aug 7 '19 at 21:26
• @wizzwizz4 The difference, then, must be that we're not connecting points, we're connecting sets of points, so we must not be deciding which points in each set to connect in an optimal way. Zoom into the text example again, the improvements you can see mostly have to do with which parts of each shape are connected. Aug 8 '19 at 18:12

# C, score 2.397x10^38

Man this took way too long to do, most likely due to my choice of language. I got the algorithm working fairly early, but ran into a lot of problems with memory allocation (couldn't recursively free stuff due to stack overflows, leak sizes were huge).

Still! It beats the other entry on every test case, and might even be optimal gets pretty close or exactly optimal solutions a lot of the time.

Anyway, here's the code:

#include <stdlib.h>
#include <stdio.h>
#include <stdbool.h>
#include <string.h>

#define WHITE 'W'
#define BLACK 'B'
#define RED   'R'

typedef struct image {
int w, h;
char* buf;
} image;

typedef struct point {
int x, y;
struct point *next;
struct point *parent;
} point;

typedef struct shape {
point* first_point;
point* last_point;

struct shape* next_shape;
} shape;

typedef struct storage {
point* points;
size_t points_size;
size_t points_index;

shape* shapes;
size_t shapes_size;
size_t shapes_index;
} storage;

char getpx(image* img, int x, int y) {
if (0>x || x>=img->w || 0>y || y>=img->h) {
return WHITE;
} else {
return img->buf[y*img->w+x];
}
}

storage* create_storage(int w, int h) {
storage* ps = (storage*)malloc(sizeof(storage));

ps->points_size = 8*w*h;
ps->points = (point*)calloc(ps->points_size, sizeof(point));
ps->points_index = 0;

ps->shapes_size = 2*w*h;
ps->shapes = (shape*)calloc(ps->shapes_size, sizeof(shape));
ps->shapes_index = 0;

return ps;
}

void free_storage(storage* ps) {
if (ps != NULL) {
if (ps->points != NULL) {
free(ps->points);
ps->points = NULL;
}
if (ps->shapes != NULL) {
free(ps->shapes);
ps->shapes = NULL;
}
free(ps);
}
}

point* alloc_point(storage* ps) {
if (ps->points_index == ps->points_size) {
printf("WHOAH THERE BUDDY SLOW DOWN\n");
/*// double the size of the buffer
point* new_buffer = (point*)malloc(ps->points_size*2*sizeof(point));
// need to change all existing pointers to point to new buffer
long long int pointer_offset = (long long int)new_buffer - (long long int)ps->points;
for (size_t i=0; i<ps->points_index; i++) {
new_buffer[i] = ps->points[i];
if (new_buffer[i].next != NULL) {
new_buffer[i].next += pointer_offset;
}
if (new_buffer[i].parent != NULL) {
new_buffer[i].parent += pointer_offset;
}
}

for(size_t i=0; i<ps->shapes_index; i++) {
if (ps->shapes[i].first_point != NULL) {
ps->shapes[i].first_point += pointer_offset;
}
if (ps->shapes[i].last_point != NULL) {
ps->shapes[i].last_point += pointer_offset;
}
}

free(ps->points);
ps->points = new_buffer;
ps->points_size = ps->points_size * 2;*/
}
point* out = &(ps->points[ps->points_index]);
ps->points_index += 1;
return out;
}

shape* alloc_shape(storage* ps) {
/*if (ps->shapes_index == ps->shapes_size) {
// double the size of the buffer
shape* new_buffer = (shape*)malloc(ps->shapes_size*2*sizeof(shape));
long long int pointer_offset = (long long int)new_buffer - (long long int)ps->shapes;
for (size_t i=0; i<ps->shapes_index; i++) {
new_buffer[i] = ps->shapes[i];
if (new_buffer[i].next_shape != NULL) {
new_buffer[i].next_shape += pointer_offset;
}
}
free(ps->shapes);
ps->shapes = new_buffer;
ps->shapes_size = ps->shapes_size * 2;
}*/
shape* out = &(ps->shapes[ps->shapes_index]);
ps->shapes_index += 1;
return out;
}

shape floodfill_shape(image* img, storage* ps, int x, int y, char* buf) {
// not using point allocator for exploration stack b/c that will overflow it

point* stack = (point*)malloc(sizeof(point));
stack->x = x;
stack->y = y;
stack->next = NULL;
stack->parent = NULL;

point* explored = NULL;
point* first_explored;
point* next_explored;

while (stack != NULL) {
int sx = stack->x;
int sy = stack->y;
stack = stack->next;

buf[sx+sy*img->w] = 1; // mark as explored

next_explored = alloc_point(ps);
next_explored->x = sx;
next_explored->y = sy;
next_explored->next = NULL;
next_explored->parent = NULL;

if (explored != NULL) {
explored->next = next_explored;
} else {
first_explored = next_explored;
}
explored = next_explored;

for (int dy=-1; dy<2; dy++) {
for (int dx=-1; dx<2; dx++) {
if (dy != 0 || dx != 0) {
int nx = sx+dx;
int ny = sy+dy;
if (getpx(img, nx, ny) == WHITE || buf[nx+ny*img->w]) {
// skip adding point to fringe
} else {
// push point to top of stack
point* new_point = (point*)malloc(sizeof(point));
new_point->x = nx;
new_point->y = ny;
new_point->next = stack;
new_point->parent = NULL;

stack = new_point;
}
}
}
}
}

/*if (getpx(img, x, y) == WHITE || buf[x+y*img->w]) {
return (shape){NULL, NULL, NULL};
} else {
buf[x+y*img->w] = 1;

shape e  = floodfill_shape(img, ps, x+1, y,   buf);
shape ne = floodfill_shape(img, ps, x+1, y+1, buf);
shape n  = floodfill_shape(img, ps, x,   y+1, buf);
shape nw = floodfill_shape(img, ps, x-1, y+1, buf);
shape w  = floodfill_shape(img, ps, x-1, y,   buf);
shape sw = floodfill_shape(img, ps, x-1, y-1, buf);
shape s  = floodfill_shape(img, ps, x,   y-1, buf);
shape se = floodfill_shape(img, ps, x+1, y-1, buf);

point *p = alloc_point(ps);
p->x = x;
p->y = y;
p->next = NULL;
p->parent = NULL;

shape o = (shape){p, p, NULL};
if (e.first_point != NULL) {
o.last_point->next = e.first_point;
o.last_point = e.last_point;
}
if (ne.first_point != NULL) {
o.last_point->next = ne.first_point;
o.last_point = ne.last_point;
}
if (n.first_point != NULL) {
o.last_point->next = n.first_point;
o.last_point = n.last_point;
}
if (nw.first_point != NULL) {
o.last_point->next = nw.first_point;
o.last_point = nw.last_point;
}
if (w.first_point != NULL) {
o.last_point->next = w.first_point;
o.last_point = w.last_point;
}
if (sw.first_point != NULL) {
o.last_point->next = sw.first_point;
o.last_point = sw.last_point;
}
if (s.first_point != NULL) {
o.last_point->next = s.first_point;
o.last_point = s.last_point;
}
if (se.first_point != NULL) {
o.last_point->next = se.first_point;
o.last_point = se.last_point;
}

return o;
}*/

shape out = {first_explored, explored, NULL};

return out;
}

shape* create_shapes(image* img, storage* ps) {
shape* first_shape = NULL;
shape* last_shape = NULL;
int num_shapes = 0;
for (int y=0; y<img->h; y++) {
for (int x=0; x<img->w; x++) {
if (getpx(img, x, y) != WHITE && !(added_buffer[x+y*img->w])) {
shape* alloced_shape = alloc_shape(ps);
*alloced_shape = floodfill_shape(img, ps, x, y, added_buffer);

if (first_shape == NULL) {
first_shape = alloced_shape;
last_shape = alloced_shape;
} else if (last_shape != NULL) {
last_shape->next_shape = alloced_shape;
last_shape = alloced_shape;
}

num_shapes++;
}
}
}

return first_shape;
}

void populate_buf(image* img, shape* s, char* buf) {
point* p = s->first_point;

while (p != NULL) {
buf[p->x+p->y*img->w] = 1;
p = p->next;
}
}

bool expand_frontier(image* img, storage* ps, shape* prev_frontier, shape* next_frontier, char* buf) {
point* p = prev_frontier->first_point;
point* n = NULL;

bool found = false;

size_t starting_points_index = ps->points_index;

while (p != NULL) {
for (int dy=-1; dy<2; dy++) {
for (int dx=-1; dx<2; dx++) {
if (dy != 0 || dx != 0) {
int nx = p->x+dx;
int ny = p->y+dy;
if ((0<=nx && nx<img->w && 0<=ny && ny<img->h) // in bounds
&& !buf[nx+ny*img->w]) {               // not searched yet
buf[nx+ny*img->w] = 1;
if (getpx(img, nx, ny) != WHITE) {
// found a new shape!
ps->points_index = starting_points_index;
n = alloc_point(ps);
n->x = nx;
n->y = ny;
n->next = NULL;
n->parent = p;
found = true;
goto __expand_frontier_fullbreak;
} else {
// need to search more
point* f = alloc_point(ps);
f->x = nx;
f->y = ny;
f->next = n;
f->parent = p;
n = f;
}
}
}
}}

p = p->next;
}
__expand_frontier_fullbreak:
p = NULL;
point* last_n = n;
while (last_n->next != NULL) {
last_n = last_n->next;
}

next_frontier->first_point = n;
next_frontier->last_point = last_n;

return found;
}

void color_from_frontier(image* img, point* frontier_point) {
point* p = frontier_point->parent;

while (p->parent != NULL) { // if everything else is right,
// a frontier point should come in a chain of at least 3
// (f point (B) -> point to color (W) -> point in shape (B) -> NULL)
img->buf[p->x+p->y*img->w] = RED;
p = p->parent;
}
}

int main(int argc, char** argv) {
if (argc < 3) {
printf("Error: first argument must be filename to load, second argument filename to save to.\n");
return 1;
}

char* fname = argv;
FILE* fp = fopen(fname, "r");

if (fp == NULL) {
printf("Error opening file \"%s\"\n", fname);
return 1;
}

int w, h;
w = 0;
h = 0;
fscanf(fp, "%d %d\n", &w, &h);

if (w==0 || h==0) {
printf("Error: invalid width/height specified\n");
return 1;
}

char* buf = (char*)malloc(sizeof(char)*w*h+1);
fgets(buf, w*h+1, fp);
fclose(fp);

image img = (image){w, h, buf};

int nshapes = 0;
storage* ps = create_storage(w, h);

while (nshapes != 1) {
// main loop, do processing step until one shape left
ps->points_index = 0;
ps->shapes_index = 0;

nshapes = 0;
while (pt != NULL) {
pt = pt->next_shape;
nshapes++;
}
if (nshapes % 1024 == 0) {
printf("shapes left: %d\n", nshapes);
}
if (nshapes == 1) {
}

shape* frontier = alloc_shape(ps);
// making a copy so we can safely free later
point* ffp = NULL;
point* flp = NULL;
while (p != NULL) {
if (ffp == NULL) {
ffp = alloc_point(ps);
ffp->x = p->x;
ffp->y = p->y;
ffp->next = NULL;
ffp->parent = NULL;
flp = ffp;
} else {
point* fnp = alloc_point(ps);
fnp->x = p->x;
fnp->y = p->y;
fnp->next = NULL;
fnp->parent = NULL;

flp->next = fnp;
flp = fnp;
}

p = p->next;
}
frontier->first_point = ffp;
frontier->last_point = flp;
frontier->next_shape = NULL;

char* visited_buf = (char*)calloc(img.w*img.h+1, sizeof(char));
populate_buf(&img, frontier, visited_buf);

shape* new_frontier = alloc_shape(ps);
new_frontier->first_point = NULL;
new_frontier->last_point = NULL;
new_frontier->next_shape = NULL;

while (!expand_frontier(&img, ps, frontier, new_frontier, visited_buf)) {
frontier->first_point = new_frontier->first_point;
frontier->last_point = new_frontier->last_point;
new_frontier->next_shape = frontier;
}

free(visited_buf);
color_from_frontier(&img, new_frontier->first_point);
img = img;
}

free_storage(ps);

char* outfname = argv;
fp = fopen(outfname, "w");

if (fp == NULL) {
printf("Error opening file \"%s\"\n", outfname);
return 1;
}

fprintf(fp, "%d %d\n", img.w, img.h);
fprintf(fp, "%s", img.buf);

free(img.buf);

fclose(fp);

return 0;
}


Tested on: Arch Linux, GCC 9.1.0, -O3

This code takes input/output in a custom file I call "cppm" (because it's like a condensed version of the classic PPM format). A python script to convert to/from it is below:

from PIL import Image

BLACK='B'
WHITE='W'
RED  ='R'

def image_to_cppm(infname, outfname):
outfile = open(outfname, 'w')
im = Image.open(infname)

w, h = im.width, im.height
outfile.write(f"{w} {h}\n")
for y in range(h):
for x in range(w):
r, g, b, *_ = im.getpixel((x, y))
if r==0 and g==0 and b==0:
outfile.write(BLACK)
elif g==0 and b==0:
outfile.write(RED)
else:
outfile.write(WHITE)
outfile.write("\n")
outfile.close()
im.close()

def cppm_to_image(infname, outfname):
infile = open(infname, 'r')

w, h = int(w), int(h)

im = Image.new('RGB', (w, h), color=(255, 255, 255))

for y in range(h):
for x in range(w):
if c==BLACK:
im.putpixel((x,y), (0, 0, 0))
elif c==RED:
im.putpixel((x,y), (255, 0, 0))

infile.close()
im.save(outfname)
im.close()

if __name__ == "__main__":
import sys
if len(sys.argv) < 3:
print("Error: must provide 2 files to convert, first is from, second is to")

infname = sys.argv
outfname = sys.argv

if not infname.endswith("cppm") and outfname.endswith("cppm"):
image_to_cppm(infname, outfname)
elif infname.endswith("cppm") and not outfname.endswith("cppm"):
cppm_to_image(infname, outfname)
else:
print("didn't do anything, exactly one file must end with .cppm")



## Algorithm explanation

How this algorithm works is that it starts by finding all the connected shapes in the image, including red pixels. It then takes the first one and expands its frontier one pixel at a time until it encounters another shape. It then colors all the pixels from the touching to the original shape (using the linkedlist it made along the way to keep track). Finally, it repeats the process, finding all the new shapes created, until there is only one shape left.

## Image gallery

Testcase 1, 183 pixels Testcase 2, 140 pixels Testcase 3, 244 pixels Testcase 4, 42 pixels Testcase 5, 622 pixels Testcase 6, 1 pixel Testcase 7, 104 pixels Testcase 8, 2286 pixels Testcase 9, 22 pixels Testcase 10, 31581 pixels Testcase 11, 21421 pixels Testcase 12, 5465 pixels Testcase 13, 4679 pixels Testcase 14, 7362 pixels • Nice work! Seems very efficient, although I can imagine a few shapes with slightly more optimal solutions: Testcase 3 (4 dots in a square), for example, I've (manually) gotten as low as 175 (a red X), not sure how I'd force that via algorithm. Aug 7 '19 at 17:09

# Python 3: 1.7x10^42 1.5x10^41

Using Pillow, numpy and scipy.

Images are assumed to be in an images folder located in the same directory as the script.

Disclaimer: It takes a long time to process all the images.

## Code

import sys
import os

from PIL import Image
import numpy as np
import scipy.ndimage

def obtain_groups(image, threshold, structuring_el):
"""
Obtain isles of unconnected pixels via a threshold on the R channel
"""
image_logical = (image[:, :, 1] < threshold).astype(np.int)
return scipy.ndimage.measurements.label(image_logical, structure=structuring_el)

def swap_colors(image, original_color, new_color):
"""
Swap all the pixels of a specific color by another color
"""
r1, g1, b1 = original_color  # RGB value to be replaced
r2, g2, b2 = new_color  # New RGB value
red, green, blue = image[:, :, 0], image[:, :, 1], image[:, :, 2]
mask = (red == r1) & (green == g1) & (blue == b1)
image[:, :, :3][mask] = [r2, g2, b2]
return image

def main(image_path=None):
images = os.listdir("images")
f = open("results.txt", "w")

if image_path is not None:
images = [image_path]

for image_name in images:
im = Image.open("images/"+image_name).convert("RGBA")
image = np.array(im)

image = swap_colors(image, (255, 255, 255), (255, 0, 0))

# create structuring element to determine unconnected groups of pixels in image
s = scipy.ndimage.morphology.generate_binary_structure(2, 2)

for i in np.ndindex(image.shape[:2]):
# skip black pixels
if sum(image[i, i]) == 255:
continue
image[i, i] = [255, 255, 255, 255]
# label the different groups, considering diagonal connections as valid
groups, num_groups = obtain_groups(image, 255, s)
if num_groups != 1:
image[i, i] = [255, 0, 0, 255]
# Show percentage
print((i + i*im.size)/(im.size*im.size))

# Number of red pixels
red_p = 0
for i in np.ndindex(image.shape[:2]):
j = (im.size - i - 1, im.size - i - 1)
# skip black and white pixels
if sum(image[j, j]) == 255 or sum(image[j, j]) == 255*4:
continue
image[j, j] = [255, 255, 255, 255]
# label the different groups, considering diagonal connections as valid
groups, num_groups = obtain_groups(image, 255, s)
if num_groups != 1:
image[j, j] = [255, 0, 0, 255]
# Show percentage
print((j + j*im.size)/(im.size*im.size))
red_p += (sum(image[j, j]) == 255*2)

print(red_p)
f.write("r_"+image_name+": "+str(red_p)+"\n")

im = Image.fromarray(image)
im.show()
im.save("r_"+image_name)
f.close()

if __name__ == "__main__":
if len(sys.argv) == 2:
main(sys.argv)
else:
main()


## Explanation

Trivial solution. We begin by changing the color of all the white pixels in an image to red. By doing this, it is guaranteed that all the elements (any isle of black pixels) are connected.

Then, we iterate over all the pixels in the image starting from the top left corner and moving right and down. For every red pixel we find we change its color to white. If after this change of color there is still only one element (an element being now any isle of black and red pixels), we leave the pixel white and move on to the next pixel. However, if after the color change from red to white the number of elements is bigger than one, we leave the pixel red and move on to the next pixel.

Update

As it can be seen (and expected) the connections obtained by only using this method show a regular pattern and in some cases, such as in the 6th and 11th images, there are unnecessary red pixels.

This extra red pixels can be easily removed by iterating again over the image and performing the same operations as explained above but from the bottom right corner to the top left corner. This second pass is much faster since the amount of red pixels that have to be checked.

## Results

The images which are modified after the second pass are listed twice to show the differences.

Number of red pixels: 18825

Number of red pixels: 334 Number of red pixels: 1352 Number of red pixels: 20214

Number of red pixels: 47268

Number of red pixels: 63 27 Number of red pixels: 17889 Number of red pixels: 259 Number of red pixels: 6746 Number of red pixels: 586

Number of red pixels: 9 1 Number of red pixels: 126 Number of red pixels: 212 Number of red pixels: 683

Score computation:

(1+6746) * (1+126) * (1+259) * (1+17889) * (1+334) * (1+586) * (1+18825) * (1+9) * (1+683) * (1+1352) * (1+20214) * (1+212) * (1+63) * (1+47268) = 1778700054505858720992088713763655500800000 ~ 1.7x10^42

Updated score computation after adding second pass:

(1+ 18825) * (1+ 1352) * (1+ 20214) * (1+ 47268) * (1+ 27) * (1+ 17889) * (1+ 6746) * (1+ 586) * (1+ 1) * (1+ 126) * (1+ 212) * (1+ 334) * (1+259) * (1+683) = 155636254769262638086807762454319856320000 ~ 1.5x10^41

• Nice work. Looks like we may need to score this one in scientific notation: 1.7x10^42 Aug 6 '19 at 17:18