I have used an [L-System](https://en.wikipedia.org/wiki/L-system) from [_The Algorithmic Beauty of Plants_](http://algorithmicbotany.org/papers/#abop) - I highly recommend it. [![ abop lsystem ][1]][1] There's a 2D L-System in the WFR [`ResourceFunction["LSystem"]`](https://resources.wolframcloud.com/FunctionRepository/resources/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] ``` [![tree][2]][2] The resulting region is filled with a large number of random points and `BinCounts` does the binning so we can get an `Image3D`: ``` ru = RegionUnion @@ tubes; pts = RandomPoint[ru, 50000]; Image3D[BinCounts[pts, {0, 800, 8}, {-200, 200, 5}, {-200, 200, 5}]] ``` [![tree image3d][3]][3] [1]: https://i.sstatic.net/kUx14l.jpg [2]: https://i.sstatic.net/aGRtW.png [3]: https://i.sstatic.net/DHlVK.png