Skip to main content

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


I'm generating some 3D models of planetary nebulae and supernova remnants for Celestia, a free OpenGL astronomy software.


Currently, I know how to do it with random points inside a spherical shell. However, I'm still unable to generate filaments like those in the Crab nebula, and I need some help on this. You can see some of the models I've created there : Nebulae models for Celestia.


Someone has a suggestion about how to generate random filaments with random points only ?


Here's a small part of the code I'm currently using to generate a shell nebula :


X[r_, u_, phi_] := Oblateness r Sqrt[1 - u^2] Cos[phi]
Y[r_, u_, phi_] := Oblateness r Sqrt[1 - u^2] Sin[phi]
Z[r_, u_] := r u

Shell[r_, u_, phi_] := {X[r, u, phi], Y[r, u, phi], Z[r, u]}


Coords[n_] := SetPrecision[Flatten[Table[Shell[#1, #2, #3] &[
Random[BetaDistribution[alpha, beta]],
Random[Real, {-1, 1}],
Random[Real, {0, 2Pi}]
], {n}], 0], CoordinatesPrecision]

This code defines "n" points inside an oblate sphere. If "n" is large, I get an uniform distribution of points, without any internal filaments-like structures.


How can I distribute the "n" points so they form "N" filaments inside the sphere, of random lenght and randomly oriented ? There should be some parameters which specify the mean number "p" of points for each filament, so approximately n = p * N.


A picture of the Crab nebula (Messier 1)



EDIT 1


Just some precision : I would like to reproduce very qualitatively the Crab nebula as a 3D object made of points, with definite cartesian coordinates X, Y, Z. The code should be compatible with Mathematica 7.0.


The features which are desired are the long filaments structures inside the nebula.


Ideally, I would like to define a statistical distribution of variables X, Y, Z that could generate some random filaments with voids between them.




EDIT 2


Here's a part of the code I'm now using to generate the models (the rest of the code isn't relevant here). Thanks a lot to all who responded, and thanks to Simon Woods, from whom that code was done ! I still have some issues, however (see below) :


InternalColor := RGBColor[0.2, 0.55, 0.8, 0.6]; (* color at the center of the nebula *)
MiddleColor := RGBColor[0.4, 0.6, 0.4, 1]; (* color of transition to the exterior part *)
ExternalColor := RGBColor[0.9, 0.3, 0.4, 0.4]; (* color of the exterior part *)

MinRadius := 0.00; (* min radius of distribution *)
MaxRadius := 1.00; (* max radius of distribution *)
MinSprite := 0.003; (* min radius of sprites *)
MaxSprite := 0.07; (* max radius of sprites *)
Oblateness := 0.8; (* oblateness of the spherical distribution *)
CoordinatesPrecision := 8;
NumberOfVoids := 1000;
NumberOfPoints := 20000;

SpriteSize[r_] := MinSprite + (MaxSprite - MinSprite)(r - MinRadius)/(MaxRadius - MinRadius);


voidpts = Select[RandomReal[{-1, 1}, {NumberOfVoids, 3}], MinRadius <= Norm[#/{1, 1, Oblateness}] <= MaxRadius &];
pts = Select[RandomReal[{-1, 1}, {NumberOfPoints, 3}], MinRadius <= Norm[#/{1, 1, Oblateness}] <= MaxRadius &];
nf = Nearest[voidpts];
DistributeDefinitions[nf];

pts = ParallelMap[Nest[0.9975 (# + 0.01 (# - First@nf[#])) &, #, 100] &, pts];

SpriteColor = Blend[{InternalColor, MiddleColor, ExternalColor}, #] &;
PlotColor = ColorData["SunsetColors"];


Graphics3D[{PointSize[0.005], {PlotColor[Norm[1.3 #]], Sphere[#, 0.005]} & /@ pts}, Boxed -> False, Background -> Black, Lighting -> "Neutral", SphericalRegion -> True]

Filaments = Join[#, {SpriteSize[Norm[#]]}, List@@SpriteColor[1.3 Norm[#]]] &/@pts;
FilamentsData := SetPrecision[Flatten[Filaments, 0], CoordinatesPrecision]

This code is slow. Is there a way to improve it ?


More importantly, I'm having an issue with the number of points generated at the intersection of several filaments : there's too much points accumulated there, and this is a problem for rendering in Celestia (too much sprites at the same location is ugly). Is there a way to reduce or dilute these points ?


Very important : I need to add a constraint on the shortest distance between two points : it shouldn't be smaller than the local sprite size. How can I add this constraint to the iteration process ?



Answer




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]

enter image description here


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}]

Comments

Popular posts from this blog

functions - Get leading series expansion term?

Given a function f[x] , I would like to have a function leadingSeries that returns just the leading term in the series around x=0 . For example: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x)] x and leadingSeries[(1/x + 2 + (1 - 1/x^3)/4)/(4 + x)] -(1/(16 x^3)) Is there such a function in Mathematica? Or maybe one can implement it efficiently? EDIT I finally went with the following implementation, based on Carl Woll 's answer: lds[ex_,x_]:=( (ex/.x->(x+O[x]^2))/.SeriesData[U_,Z_,L_List,Mi_,Ma_,De_]:>SeriesData[U,Z,{L[[1]]},Mi,Mi+1,De]//Quiet//Normal) The advantage is, that this one also properly works with functions whose leading term is a constant: lds[Exp[x],x] 1 Answer Update 1 Updated to eliminate SeriesData and to not return additional terms Perhaps you could use: leadingSeries[expr_, x_] := Normal[expr /. x->(x+O[x]^2) /. a_List :> Take[a, 1]] Then for your examples: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x), x] leadingSeries[Exp[x], x] leadingSeries[(1/x + 2 + (1 - 1/x...

How to thread a list

I have data in format data = {{a1, a2}, {b1, b2}, {c1, c2}, {d1, d2}} Tableform: I want to thread it to : tdata = {{{a1, b1}, {a2, b2}}, {{a1, c1}, {a2, c2}}, {{a1, d1}, {a2, d2}}} Tableform: And I would like to do better then pseudofunction[n_] := Transpose[{data2[[1]], data2[[n]]}]; SetAttributes[pseudofunction, Listable]; Range[2, 4] // pseudofunction Here is my benchmark data, where data3 is normal sample of real data. data3 = Drop[ExcelWorkBook[[Column1 ;; Column4]], None, 1]; data2 = {a #, b #, c #, d #} & /@ Range[1, 10^5]; data = RandomReal[{0, 1}, {10^6, 4}]; Here is my benchmark code kptnw[list_] := Transpose[{Table[First@#, {Length@# - 1}], Rest@#}, {3, 1, 2}] &@list kptnw2[list_] := Transpose[{ConstantArray[First@#, Length@# - 1], Rest@#}, {3, 1, 2}] &@list OleksandrR[list_] := Flatten[Outer[List, List@First[list], Rest[list], 1], {{2}, {1, 4}}] paradox2[list_] := Partition[Riffle[list[[1]], #], 2] & /@ Drop[list, 1] RM[list_] := FoldList[Transpose[{First@li...

front end - keyboard shortcut to invoke Insert new matrix

I frequently need to type in some matrices, and the menu command Insert > Table/Matrix > New... allows matrices with lines drawn between columns and rows, which is very helpful. I would like to make a keyboard shortcut for it, but cannot find the relevant frontend token command (4209405) for it. Since the FullForm[] and InputForm[] of matrices with lines drawn between rows and columns is the same as those without lines, it's hard to do this via 3rd party system-wide text expanders (e.g. autohotkey or atext on mac). How does one assign a keyboard shortcut for the menu item Insert > Table/Matrix > New... , preferably using only mathematica? Thanks! Answer In the MenuSetup.tr (for linux located in the $InstallationDirectory/SystemFiles/FrontEnd/TextResources/X/ directory), I changed the line MenuItem["&New...", "CreateGridBoxDialog"] to read MenuItem["&New...", "CreateGridBoxDialog", MenuKey["m", Modifiers-...