Here's the code. The cover image used m=20, and the final ImageSize was set to 5200 (the highest value I could get without risking a crash). I exported the image as a PNG, and used Photoshop's shift-tilt filter to introduce the blurred edges. I could have done that in mma, of course, but it's hard to argue with the simplicity of the photoshop filter. If you run this with something small, like m=4, you can easily grab and rotate the resulting image to see the pyramid structure.
Clear[s, e, n, w, dir, tab, m, \[Delta], positions, allPrimes,
twinPrimes, sph, ulam];
m = 4; (* m is the no. of spiral "laps", so (2m)^2 numbers will be \
displayed *)
\[Delta] = .1; (* \[Delta] = vertical depth added to \
pyramid for each number *)
dir = {s, e, n,
w}; (* directions: south, east, north, west *)
tab = {};
Do[dir = RotateLeft[dir]; AppendTo[tab, Table[First[dir], {k}]];
dir = RotateLeft[dir];
AppendTo[tab, Table[First[dir], {k}]], {k, 1, 2 m}];
positions =
FoldList[Plus, {0, 0, 0},
Take[Flatten[tab], 4 m^2 - 1] /. {e -> {\[Delta], -1, 0},
n -> {\[Delta], 0, 1}, w -> {\[Delta], 1, 0},
s -> {\[Delta], 0, -1}}];
allPrimes = PrimeQ[Range[4 m^2]];
twinPrimes =
MapThread[
Or, {MapThread[
And, {allPrimes, Join[{False, False}, Drop[allPrimes, -2]]}],
MapThread[
And, {allPrimes, Join[Drop[allPrimes, 2], {False, False}]}]}];
sph[k_, pos_, highlight_] :=
Translate[
SphericalPlot3D[.3, {\[Theta], 0, Pi}, {\[Phi], 0, 2 Pi},
Mesh -> None,
PlotStyle ->
Texture[Graphics[
Style[Text[k, {0, 0}, {0, 0}, {0, 1}], 128,
FontFamily -> "Optima"], AspectRatio -> 2,
Background ->
If[highlight, RGBColor[1.`, 1.`, 0.`], RGBColor[.5`, .5`, .5`]],
ImageSize -> 400]],
TextureCoordinateFunction -> ({#4, #5} &)][[1]], pos];
ulam = Graphics3D[{
{Gray, Tube[positions, .05]},
MapThread[sph, {Range[4 m^2], positions, twinPrimes}]
}, ViewPoint -> {-1.8, 0, 0}, PlotRangePadding -> None,
Boxed -> False, Lighting -> "Neutral", ImageSize -> 600,
Background -> Black]
