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