# How to generate capillaries diagram?

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?

• You may want to look at this visualization.
– Syed
Commented Dec 18, 2021 at 15:53
• You can not create a 3D picture from a 2D picture. The former has a lot more information than the latter. Commented Dec 19, 2021 at 9:13
• @Syed Thanks for the suggestion!
– SAC
Commented Dec 21, 2021 at 2:45
• @DanielHuber thanks for the answer.
– SAC
Commented Dec 21, 2021 at 2:45
• Related: 3d-tree-in-mathematica Commented Dec 23, 2021 at 3:23

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


• Amazing! Congratulations and thank you so much
– SAC
Commented Dec 23, 2021 at 9:58