Code:

import math,sympy
Size = 1025 # must be an odd number
offset = Size>>1 # so the center is at (512,512) if Size==1025
def get_x_y(n): # this function is from http://oeis.org/A235913
sr = int(math.sqrt(n-1)) # OK for relatively small n's
sr = sr-1+(sr&1)
rm = n-sr*sr
d = (sr+1)/2
if rm<=sr+1:
return -d+rm, d
if rm<=sr*2+2:
return d, d-(rm-(sr+1))
if rm<=sr*3+3:
return d-(rm-(sr*2+2)), -d
return -d, -d+rm-(sr*3+3)
outfile = open("outfile_42.ppm", "wb")
outfile.write("P6\n"+str(Size)+" "+str(Size)+"\n255\n")
data = [4] * (Size*Size)
p = runLength = 1
print 'Starting the main loop...'
for i in range(Size*Size): # runs ~2 seconds if Size==1025
runLength -= 1
if runLength==0:
p = sympy.nextprime(p+1)
colorCode = (p & 7) >> 1 # extract 3 bits, discard the last bit, it's always a 1
runLength = (p >> 3) % 42 + 1 # 42 is a magic number
x,y = get_x_y(i+1)
data[x+offset+(y+offset)*Size] = colorCode
print 'Writing image file...'
for x in range(Size): # also ~2 seconds if Size==1025
for y in range(Size):
rgb = [0,0,0]
v = data[x+y*Size]
if v<4:
numNeighbors = 0
for i in range(x-3, x+3 +1):
for j in range(y-3, y+3 +1):
if i>=0 and j>=0 and i<Size and j<Size and data[i+j*Size]<4: numNeighbors+=1
if numNeighbors > 1: # apply some noise reduction. > 0 to skip noise reduction
if v<3: rgb[v] = 255
elif v==3: rgb[0] = rgb[1] = rgb[2] = 255
for c in rgb: outfile.write(chr(c))
outfile.close()