How to generate capillaries diagram?

Question

I would like to generate a 3D matrix like the 3D figure below.

How to generate a tubular shell branching process in 3D like that image? (veins dynamic from cellular automata?)

n=128; (*dimension of matrix*)
f[i_,j_,k_]:= (*rule to generate the branches tube*)
s = SparseArray[{{i_, j_,k_} -> f[i, j, k]}, {n, n, n}]; (*matrix*)
t = Tube[{s}, {0.6, 0.4}]; (*tube*)

Could anybody help me to obtain the matrix, please?

Answer 1

I have used an L-System from The Algorithmic Beauty of Plants - I highly recommend it.

There's a 2D L-System in the WFR ResourceFunction["LSystem"], but unfortunately not 3D yet, so I made one for this answer. This L-System is parametric so I'm using RegularExpression to extract the parameters in ( ... ) parentheses and update them during replacement.

proc[s_, params_] := ToString@ToExpression[s]
lsys[rules_, axiom_, n_] :=
 Nest[StringReplace[#, rules] &, axiom, n]

convert[str_] := StringReplace[StringCases[str,
   RegularExpression[".\\(.+?\\)"] |
    "[" | "]" | RegularExpression["."]
   ], {"(" -> "[", ")" -> "]",
   "[" -> "branch[]", "]" -> "complete[]",
   "F" -> "forward",
   "/" -> "rollright",
   "\\" -> "rollleft",
   "&" -> "pitchdown",
   "∧" -> "pitchup",
   "+" -> "turnleft",
   "-" -> "turnright",
   "$" -> "vertical[]",
   "!" -> "diameter"
   }]

(* parameters *)
θ = 94.74;
ϕ = 132.63;
a = 18.95;
e = 1.109;
v = 1.7;
axiom = "!(1)F(200.)/(45.)A";
params = #[[1]] -> ToString[#[[2]]] & /@ {"θ" -> θ, "ϕ" -> ϕ, "a" -> a, "e" -> e, "v" -> v};
rules = {
   "A" -> StringReplace["!(v)F(50)[&(a)F(50)A]/(θ)[&(a)F(50)A]/(ϕ)[&(a)F(50)A]", params],
   RegularExpression@"F\\((.+?)\\)" :> "F(" <> proc["e*$1", params] <> ")",
   RegularExpression@"!\\((.+?)\\)" :> "!(" <> proc["v*$1", params] <> ")"
};

stack = CreateDataStructure["Stack"];
branches = CreateDataStructure["LinkedList"];
state = <|"frame" -> IdentityMatrix[3], "position" -> {0, 0, 0}, "diameter" -> 10|>;
branch[] := stack["Push", state];
complete[] := If[! stack["EmptyQ"], state = stack["Pop"]];
forward[x_] := 
  With[{prev = state["position"]}, 
   state["position"] = state["position"] + state["frame"][[1]]*x;
   branches["Append", Cylinder[{prev, state["position"]}, state["diameter"]]]];
diameter[x_] := state["diameter"] = x;
turnleft[x_] := state["frame"] = RotationMatrix[x °, {0, 0, 1}].state["frame"];
turnright[x_] := state["frame"] = RotationMatrix[-x °, {0, 0, 1}].state["frame"];
pitchdown[x_] := state["frame"] = RotationMatrix[x °, {0, 1, 0}].state["frame"];
pitchup[x_] := state["frame"] = RotationMatrix[-x °, {0, 1, 0}].state["frame"];
rollleft[x_] := state["frame"] = RotationMatrix[x °, {1, 0, 0}].state["frame"];
rollright[x_] := state["frame"] = RotationMatrix[-x °, {1, 0, 0}].state["frame"];
vertical[] := state["frame"] = IdentityMatrix[3];

ToExpression /@ convert[lsys[rules, axiom, 5]];
tubes = Normal@branches;
gr = Graphics3D[tubes, ViewVertical -> {1, 0, 0}, Boxed -> False]

The resulting region is filled with a large number of random points and BinCounts does the binning and creates the matrix we can view as an Image3D:

ru = RegionUnion @@ tubes;
pts = RandomPoint[ru, 50000];
Image3D[BinCounts[pts, {0, 800, 8}, {-200, 200, 5}, {-200, 200, 5}]]