Fast surface area counting (optimized)

import collections
import functools

lines = [l.strip() for l in open('day18.txt', 'r') if l.strip()]
P = set(map(eval, lines))

def adj(p):
    for axis in range(3):
        for d in (-1, 1):
            q = list(p)
            q[axis] += d
            yield tuple(q)

P_adj = {q for p in P for q in adj(p) if q not in P}

lines = collections.defaultdict(lambda: [])
indices = {}
for p in P_adj:
    for axis in range(3):
        lines[(axis, p[:axis] + p[axis + 1:])].append(p)
for (axis, _), l in lines.items():
    l.sort(key=lambda p: p[axis])
    for i, p in enumerate(l):
        indices[(axis, p)] = i


def get_line(p, axis):
    l = lines[(axis, p[:axis] + p[axis + 1:])]
    i = indices[(axis, p)]
    return l, i

def neighbors(p):
    #Try to jump across empty space
    for axis in range(3):
        l, i = get_line(p, axis)
        if i == 0 or i == len(l) - 1:
            continue
        for d in (-1, 1):
            if not (0 <= i + d < len(l)):
                continue
            q = list(p)
            q[axis] += d
            q = tuple(q)
            if q in P:
                continue
            yield l[i + d]
    #Crawl around the boundary
    for q in adj(p):
        if q in P:
            continue
        if q in P_adj:
            yield q
        for r in adj(q):
            if r in P_adj:
                yield r

external = set()

#Find points which directly face the outside
for l in lines.values():
    external.add(l[0])
    external.add(l[-1])

#Flood fill
stack = list(external)
while stack:
    p = stack.pop()
    for q in neighbors(p):
        if q not in external:
            external.add(q)
            stack.append(q)

print(sum(1 for p in external for q in adj(p) if q in P))