How to represent flux in a shape

Question

I used Tim Laska's solution here for raytracing, and modelled 10000 photons impinging onto a cone.

shape = Cone[{{7.5, 7.5, 1}, {7.5, 7.5, 6}}, 5]

(* 10 000 photons *) 
data2 = Import[
    "https://raw.githubusercontent.com/Tb8854/overflowquestion/master/\
test", "Data"][[1]];
data2 = ToExpression[data2];

(* plot first 1000 *) 
data = data2[[1 ;; 1000]];
data = Flatten[data[[#]], 2] & /@ Range[Length[data]];
lines = Line[#] & /@ data;
Graphics3D[{lines, {Opacity[0.4], shape}}, 
 PlotRange -> {{2, 18}, {2, 18}, {-3, 8}}, Axes -> True, 
 Boxed -> False]

But, I want to see the flux, i.e. the number of rays (i.e. lines) per unit area in the material. Maybe something like this:

Where, if I a higher density of lines the cone, the redder the cone should be. I don't think the rays in the first image effectively show the spreading out of the rays when they enter the cone.

---

A Poor Solution

I've developed a poor 2D solution, beware - it takes a while (2 minutes). What it does is make slices in the cone, and check if there is an intersection. It works okay.

datastore = ConstantArray[{}, (5)];
points = ConstantArray[{}, {(5), Length[data]}];
(* prevent hitting the edge *) 
positionerror = 0.01;

For[j = 5, j >= 1, j = j - 1,
 Print[j]; 
 number = 0; 
 plane = InfinitePlane[{10, 10, 
    j + positionerror}, {{0, 1, 0}, {1, 0, 0}}];

 Monitor[For[i = 1, i <= Length[data], i++, 
   line = Line[data[[i]]];
   intersection = RegionIntersection[plane, line, shape];
   points[[j, i]] = intersection;
   If[intersection === EmptyRegion[3], 
    Continue[];, 
    number = number + 1;
    ];

   area = Area[Region[RegionIntersection[plane, shape]]];
   datastore[[j]] = {number, area, number/area};

   ], ProgressIndicator[i, {1, Length[data]}]]]


For[i = 2.01, i <= 6, i++, 
  pointsplot = 
   Select[Flatten[points[[2 ;;, ;; , 1]], 2], Length[#] == 3 &];
  pointsplot = Select[pointsplot, #[[3]] == i &];
  pointsplot = pointsplot[[;; , {1, 2}]];

  data = BinCounts[pointsplot, {2, 13, 0.1}, {2, 13, 0.1}];
  data = GaussianFilter[data, 2];
  Print[ListDensityPlot[data, ColorFunction -> "BlueGreenYellow", 
    PlotRange -> All]]
  ];

Graphics3D[{points}]

But does anyone have any ideas how to maybe represent the flux in 3D?

Answer 1

Maybe SmoothDensityHistogram could help smooth out and visualize the flux.

poincarePlane[p_] := 
 Select[Flatten[points[[2 ;;, ;; , 1]], 2], 
   Length[#] == 3 && (Abs[#[[3]] - p] < 0.01) &][[All, {1, 2}]]
Column[SmoothDensityHistogram[poincarePlane[#], 
    ColorFunction -> 
     Function[{z}, ColorData["BlueGreenYellow"][10 z]], 
    ColorFunctionScaling -> False, 
    PlotRange -> RegionBounds[shape][[1 ;; 2]]] & /@ 
  Reverse@{2.01, 3.01, 4.01, 5.01}]

You could try constructing a 3D model using Image3D as describe by C.E. in this broccoli reconstruction MSE answer here.

Update

I have not tried to optimize Tomi's code, but I increased the number of planes from 5->20.

datastore = ConstantArray[{}, (20)];
points = ConstantArray[{}, {(20), Length[data]}];
(*prevent hitting the edge*)
positionerror = 0.01;

For[j = 20, j >= 1, j = j - 1, Print[j];
 number = 0;
 plane = InfinitePlane[{10, 10, 
    j/4 + positionerror}, {{0, 1, 0}, {1, 0, 0}}];
 Monitor[For[i = 1, i <= Length[data], i++, line = Line[data[[i]]];
   intersection = RegionIntersection[plane, line, shape];
   points[[j, i]] = intersection;
   If[intersection === EmptyRegion[3], Continue[];, 
    number = number + 1;];
   area = Area[Region[RegionIntersection[plane, shape]]];
   datastore[[j]] = {number, area, number/area};], 
  ProgressIndicator[i, {1, Length[data]}]]]

Here are the images at finer z-spacing:

cuts = Table[i + 0.01, {i, 1, 5., 0.25}];
poincarePlane[p_] := 
 Select[Flatten[points[[2 ;;, ;; , 1]], 2], 
   Length[#] == 3 && (Abs[#[[3]] - p] < 0.01) &][[All, {1, 2}]]
imgs = Rasterize@
    SmoothDensityHistogram[poincarePlane[#], 
     ColorFunction -> 
      Function[{z}, ColorData["BlueGreenYellow"][20 z]], 
     ColorFunctionScaling -> False, Background -> Black, 
     PlotRange -> RegionBounds[shape][[1 ;; 2]], ImageSize -> Small, 
     FrameTicks -> None, ImagePadding -> None] & /@ Reverse@cuts

3D Reconstruction Using Image3D

Using @C.E.'s approach linked above:

imgs = MapThread[SetAlphaChannel, {imgs, Binarize /@ imgs}];
Image3D[imgs, Background -> Black, BoxRatios -> {1, 1, 1}]

Dynamic Cut Plane Visualization

The Texture shows how to create an efficient dynamic clip planes visualization.

data = Developer`ToPackedArray[Map[ImageData, imgs]];
Manipulate[
 Graphics3D[{Opacity[Dynamic[o]], Texture[data], EdgeForm[None],
   Dynamic[{Polygon[{{x, 0, 0}, {x, 1, 0}, {x, 1, 1}, {x, 0, 1}}, 
      VertexTextureCoordinates -> {{x, 0, 0}, {x, 1, 0}, {x, 1, 
         1}, {x, 0, 1}}], 
     Polygon[{{0, y, 0}, {1, y, 0}, {1, y, 1}, {0, y, 1}}, 
      VertexTextureCoordinates -> {{0, y, 0}, {1, y, 0}, {1, y, 
         1}, {0, y, 1}}],
     Polygon[{{0, 0, 1 z}, {1, 0, 1 z}, {1, 1, 1 z}, {0, 1, 1 z}}, 
      VertexTextureCoordinates -> {{0, 0, 1 z}, {1, 0, 1 z}, {1, 1, 
         1 z}, {0, 1, 1 z}}]}]}, Background -> Black, 
  RotationAction -> "Clip"], {{x, 0.5}, 0, 1}, {{y, 0.5}, 0, 
  1}, {{z, 0.5}, 0, 1}, {{o, 0.75, "opacity"}, 0, 1}, 
 ControlPlacement -> Top]

Update to Include Volume Rendering Customization

One may desire finer control over the Opacity to obtain the desired volume rendering. @Jason B's MSE answer here shows a way to achieve this using ListDensityPlot3D.

In this case, I ran @Tomi's work flow on the full dataset on 40 planes overnight and applied the following workflow to generate an Image3D.

cuts = Table[i + 0.01, {i, 1, 5., 0.125}];
pts = Transpose[
   Select[Flatten[points[[2 ;;, ;; , 1]], 2], Length[#] == 3 &]];
poincarePlane[p_] := 
 Select[Flatten[points[[2 ;;, ;; , 1]], 2], 
   Length[#] == 3 && (Abs[#[[3]] - p] < 0.01) &][[All, {1, 2}]]
imgs = Rasterize@
     SmoothDensityHistogram[poincarePlane[#], 
      ColorFunction -> Function[{z}, ColorData["GrayTones"][40 z]], 
      PlotRange -> {{1.5, 13.5}, {1.5, 13.5}}, 
      RegionFunction -> Function[{x, y, z}, z >= 0.001], 
      ColorFunctionScaling -> False, Background -> Black, 
      ImageSize -> Medium, FrameTicks -> None] & /@ Reverse@cuts;
imgs = ColorConvert[ImageCrop[#, 590] & /@ imgs, "Grayscale"];
i3d = Image3D[imgs, Background -> Black, BoxRatios -> {1, 1, 1}, 
  ColorFunction -> "XRay"]

Now, we can apply @Jason B's workflow to explore different OpacityFunctions.

list = ImageData /@ Image3DSlices[i3d];
Plot[Evaluate[(Exp[# f] - 1)/(E^# - 1) & /@ {1, 4, 8, 12}], {f, 0, 1}]
ListDensityPlot3D[list, 
   OpacityFunction -> Function[f, (Exp[# f] - 1)/(E^# - 1)], 
   ColorFunction -> ColorData["BlueGreenYellow"], 
   RegionFunction -> 
    Function[{x, y, z, f}, ! ((x >= 7.5) && (y <= 7.5))], 
   ImageSize -> 300, Background -> Black, 
   DataRange -> MinMax /@ pts] & /@ {1, 4, 8, 12}