Distribution of random points in 3D space to simulate the Crab Nebula

Here's an attempt in which I start with a set of "void points", which will be the centres of the gaps between filaments. The stars are then created as an initially random distribution, and are repeatedly nudged away from their nearest void point. Or, to look at it another way, they are attracted towards the edges of the Voronoi cells defined by the void points. There is also a contraction at each iteration to prevent stars from escaping towards infinity.

Update: My original code used ParallelMap which turned out to be considerably slower than plain old Map. Thanks to @halmir for pointing that out. I should have checked. Following that observation, a couple of other optimisations emerged. The updated code below includes simplification of the filamentation function using Expand, and performs each iteration step on all the points at once (where the original code did the full set of iterations on each point in turn). Finally, I have used ParallelTable to distribute the calculation across all the CPU cores (and this time I have checked that it is faster that way...) The filementation code now runs in a couple of seconds on my 4-core machine.

A bit more explanation as requested:

In the code below the first two lines create random points in the sphere (actually it would be better to use the OP's Coords function in the question for this bit, as it provides better control of the point density).

The function f embodies the the filamentation. At each step, each point in pts moves 1% further away from its nearest void point (the 0.01) and then 0.25% closer to the origin (the 0.9975). This is repeated 200 times for each point.

The filamentation process reduces the size of the point cloud from the initial radius of 1 to something around 0.75. In the Graphics3D output the points are coloured according to their distance from the origin, the factor of 1.3 in the colour function is simply to there to compensate for the contraction.

voidpts = Select[RandomReal[{-1, 1}, {1000, 3}], Norm[#] <= 1 &];
pts = Select[RandomReal[{-1, 1}, {10000, 3}], Norm[#] <= 1 &];
nf = Nearest[voidpts];

f = Evaluate[Expand[0.9975 (#1 + 0.01 (#1 - #2))]] &;

DistributeDefinitions[nf, f];
pts = Developer`ToPackedArray @ Drop[pts, Length[pts] ~ Mod ~ $KernelCount];

pts = Join @@ ParallelTable[Nest[f[#, (nf /@ #)[[All, 1]]] &, p, 200],
    {p, Partition[pts, Length[pts]/$KernelCount]}];

cf = ColorData["SunsetColors"];
Graphics3D[{PointSize[0.005], {cf[Norm[1.3 #]], Sphere[#, 0.005]} & /@ pts}, 
  Boxed -> False, Background -> Black, Lighting -> "Neutral", 
  SphericalRegion -> True, ViewPoint -> {0.75, -0.75, 0.75}, ViewAngle -> 0.7]

In principle you could muck around with the initial distribution of void points to get larger voids and better defined filaments in the centre, as appears to be the case in the original image. Unfortunately the code is rather slow so it's not much fun to experiment with.

Here's the gratuitous animation:

To create an ouput list of the form {{X1, Y1, Z1, R1, G1, B1}, {X2, Y2, Z2, R2, G2, B2}, {X3, Y3, Z3, R3, G3, B3}, ...} you could do something like:

coords = Join[#, List @@ cf[1.3 Norm[#]]] & /@ pts;

Update 2021

Astronomers have actually done the job for you with the real crab nebula:

http://www.cfht.hawaii.edu/en/news/HeartofTheCrab/out-outreach.mp4

This movie (be a bit patient) allows you to rotate around the observed nebulae as described here

---

If you use this code take the hessian of it and plot the map of the largest eigenvalues you get a nice filamentary map like the bottom right panel.

see this reference (specifically pp 28 of the phd).

In mathematica it can be coded as follows

 nn = 256;u = GaussianRandomField[nn, 2, Function[k, k^-4]]//GaussianFilter[#, 4] & // Chop;
  Clear[f]; f[x_, y_] = ListInterpolation[u][x, y];
  Clear[d2f]; d2f[x_, y_] = {{D[f[x, y], {x, 2}],
        D[f[x, y], x, y]}, {D[f[x, y], x, y], D[f[x, y], {y, 2}]}};
  h = Table[d2f[x, y], {x, nn}, {y, nn}];
  ee = Table[Eigenvalues[h[[i, j, All, All]]], {i, nn}, {j, nn}];
  ee = Map[Max, Abs[ee], {2}];
  ee[[1 ;; nn - 1, 1 ;; nn - 1]] // Image // ImageAdjust

producing

Removing the Abs for nn=512 yields

Varying the power law and the smoothing would allow you to produce e.g.

and if you mix the result of 3 such images

Transpose[{ee1, ee2, ee3}, {3, 1, 2}] // Image // ImageAdjust

The same applies in 3D

nn = 64; u = GaussianRandomField[nn, 3, Function[k, k^-3]]//GaussianFilter[#, 6] & // Chop; 
Clear[f]; f[x_, y_, z_] = ListInterpolation[u][x, y, z];
Clear[d2f]; d2f[x_, y_, z_] = {
  {D[f[x, y, z], x, x], D[f[x, y, z], x, y], 
   D[f[x, y, z], x, z]}, {D[f[x, y, z], x, y], D[f[x, y, z], y, y], 
   D[f[x, y, z], y, z]},
  {D[f[x, y, z], x, z], D[f[x, y, z], y, z], D[f[x, y, z], z, z]}};
h = Table[d2f[x, y, z], {x, nn}, {y, nn}, {z, nn}];

ee = Table[
   Eigenvalues[h[[i, j, k, All, All]]], {i, nn}, {j, nn}, {k, nn}];
ee = Map[Max, ee, {3}];

If we look at a slice in 3D

 ImageAdjust@Image3D[Exp[ee/StandardDeviation[Flatten[ee]]], Background -> Black]

EDIT

Let us explore another venue, just because this is how cosmologist generate initial conditions for simulations. Let us generate a displacement field corresponding to the gradient of a Gaussian Random field:

nn = 256; u = 
 GaussianRandomField[nn, 2, Function[k, k^-4]] // 
   GaussianFilter[#, 4] & // Chop;
Clear[f]; f[x_, y_] = ListInterpolation[u][x, y];
df[x_, y_] = {D[f[x, y], x], D[f[x, y], y]};
g = Table[df[x, y], {x, nn}, {y, nn}]; g /= Max[g];
grid = Outer[{#1, #2} &, Range[nn], Range[nn]];
Map[Point, grid + g*15, {2}] // 
 Graphics[{AbsolutePointSize[0.1], #}] &

and in 3D

nn = 32; u = GaussianRandomField[nn, 3, Function[k, k^-4]] //GaussianFilter[#, 4] & // Chop;
Clear[f]; f[x_, y_, z_] = ListInterpolation[u][x, y, z];
df[x_, y_, z_] = {D[f[x, y, z], x], D[f[x, y, z], y], D[f[x, y, z], z]};
g = Table[df[x, y, z], {x, nn}, {y, nn}, {z, nn}]; g /= Max[g];
grid = Table[{i, j, k}, {i, nn}, {j, nn}, {k, nn}];
Map[Point, grid + g*8, {2}] // Graphics3D[{AbsolutePointSize[0.1], #}] &

Or to make it look good borrowing Simon's graphics primitive

cf = ColorData["SunsetColors"]
Graphics3D[{PointSize[0.1], {cf[#[[2]]/nn], Sphere[#, 0.2]} & /@ pts}, 
Boxed -> False, Background -> Black, Lighting -> "Neutral", SphericalRegion -> True]

Now the above code is fast (and could be made faster while using FFT to compute the gradients).

Note that we can remove the regular low density region points as follows

 pts = MapThread[If[Norm[#1] > 1/3, 15 #1 + #2, {}] &, {g, grid}, 2] //
  Flatten[#, 1] & // Union // Rest;
 Map[Point, pts] // Graphics[{AbsolutePointSize[0.1], #}] &

EDIT 2

If I steal the set of points from Simon above (I hope he doesn't mind!) and trace its 3D skeleton I get this (So the credit to this method lies with Simon Woods and Thierry Sousbie)

or changing the so called level of persistence and increasing the number of drawn points in Simon's code:

And its really 3D :-)

Note that the original 2D image can be analysed directly by the skeleton, channel per channel,

which demonstrates that the filamentary structure is colour dependent, but I guess I am getting over-carried! :-)

COMMENT

To answer the OP, the skeleton is not currently implemented in mathematica. The closest current implementation is watershading:

nn = 512; u0 =  GaussianRandomField[nn, 2, Function[k, k^-3]]//GaussianFilter[#, 3] & // Chop;
u = u0 // GaussianFilter[#, 8] & // Chop;
Clear[f]; f[x_, y_] = ListInterpolation[u][x, y];
df[x_, y_] = {D[f[x, y], x], D[f[x, y], y]};
g = Table[df[x, y] // Sqrt[#.#] &, {x, nn}, {y, nn}];
im0 = u0 // Image;im1 = g // Image // WatershedComponents // Image;
ImageMultiply[im1, im0] // ImageAdjust

May be in version 10 ;-)

Here is the modified version of Simon Woods answer. In my machine with 2 core, ParallelMap version is slower and you could gain a little bit of speed up by using compiled version of iteration:

voidpts = Select[RandomReal[{-1, 1}, {1000, 3}], Norm[#] <= 1 &];
pts = Select[RandomReal[{-1, 1}, {10000, 3}], Norm[#] <= 1 &];
nf = Nearest[voidpts];
DistributeDefinitions[nf];

cNest = Compile[{{x, _Real, 1}, {n, _Integer}}, 
         Block[{pt = x}, 
             Do[pt = 0.9975 (pt + 0.01 (pt - First@nf[pt])), {i, n}]; 
             pt], {{nf[_], _Real, 2}}];

The following is the timing I got:

In[75]:= pt1 = 
    ParallelMap[Nest[0.9975 (# + 0.01 (# - First@nf[#])) &, #, 100] &, 
     pts]; // AbsoluteTiming
Out[75]= {98.073226, Null}

In[85]:= pt2 = 
   Map[Nest[0.9975 (# + 0.01 (# - First@nf[#])) &, #, 100] &, 
    pts]; // AbsoluteTiming
Out[85]= {7.226970, Null}

In[86]:= pt3 = Map[cNest[#, 100] &, pts]; // AbsoluteTiming
Out[86]= {5.901947, Null}

In[80]:= pt1 == pt2 == pt3
Out[80]= True

I don't see the any reason to use Sphere, so I replace with Point to gain speed:

cf = ColorData["SunsetColors"]; 
Graphics3D[{AbsolutePointSize[2], 
     Point[pt1, VertexColors -> (cf[Norm[1.3 #]] & /@ pt1)]}, 
     Boxed -> False, SphericalRegion -> True, 
     ViewPoint -> {0.75, -0.75, 0.75}, ViewAngle -> 0.7, 
     Background -> Black]

To check how each iteration goes (for fun):

ptlist = Transpose[
      Map[NestList[0.9975 (# + 0.01 (# - First@nf[#])) &, #, 200] &, 
     pts]];

Manipulate[
  Graphics3D[{AbsolutePointSize[2], 
    Point[ptlist[[i]], 
    VertexColors -> (cf[Norm[1.3 #]] & /@ ptlist[[i]])]}, 
    Boxed -> False, SphericalRegion -> True, 
    ViewPoint -> {0.75, -0.75, 0.75}, ViewAngle -> 0.7, 
    Background -> Black, PlotRange -> {-1.5, 1.5}], 
 {i, 1, Length[ptlist], 1}]