3D tree in Mathematica?

Question

Searching the web for information about the affine transformation, I found the one page, which called my attention for the tree that show and is this

but unfortunately do not give information about the algorithm to create it, I would like to ask help to make one the same or very similar, maybe someone knows where to get more information about it. Thanks in advance, here is the link to the page mentioned above

Answer 1

First, an idomatic, but slow version.

s1 = 1/GoldenRatio // N;
s2 = 1/GoldenRatio // N;
stem = {0., 0., 1.};
thickness = 0.15;
branches = Table[RotationMatrix[2. k Pi/3., {0, 0, 1}].{Cos[Pi/4.], 0., Sin[Pi/4.]}, {k, 0, 2}];
data0 = {Join[{{0., 0., 0.}}, {stem}, branches, {{thickness, 1., 0.}}]};
iteration[data_] :=
 Block[{U},
  Flatten[Table[
    U = data[[j]];
    Table[
     Join[{U[[1]] + U[[2]]}, {U[[i]]}, 
      s1 U[[3 ;; 5]].RotationMatrix[{U[[i]], U[[2]]}], {s2 U[[6]]}],
     {i, 3, 5}],
    {j, 1, Length[data]}
    ],
   1
   ]
  ]

This generates the tree structure.

result = NestList[iteration, data0, 6]; // AbsoluteTiming

(* {0.211536, Null} *)

This generates the tree plot.

t = 0.5;
colfun[x_] := ColorData["Rainbow"][t + (1 - t) x];
plot[U_] := {colfun[U[[6, 2]]],Table[Sphere[U[[1]] + t U[[2]], U[[6, 1]] (1 - t) + t s2 U[[6, 1]]], {t, 0.0, 0.9, 0.1}]};
Graphics3D[
 Flatten[plot /@ Flatten[result, 1]],
 Lighting -> "Neutral",
 Background -> Black,
 Boxed -> False,
 SphericalRegion -> True
 ]

A faster version is generated with Compile and some handcraft:

citeration2 = 
  With[{scale1 = s1, scale2 = s2, part = Compile`GetElement},
   Compile[{{U, _Real, 2}}, Block[{A, u, v, w},
     v = {part[U, 2, 1], part[U, 2, 2], part[U, 2, 3]}/Sqrt[part[U, 2, 1]^2 + part[U, 2, 2]^2 + part[U, 2, 3]^2];
     Table[
      u = {part[U, i, 1], part[U, i, 2], part[U, i, 3]}/Sqrt[part[U, i, 1]^2 + part[U, i, 2]^2 + part[U, i, 3]^2];
      w = {
        -(part[u, 3] part[v, 2]) + part[u, 2] part[v, 3],
        part[u, 3] part[v, 1] - part[u, 1] part[v, 3],
        -(part[u, 2] part[v, 1]) + part[u, 1] part[v, 2]
        };
      w = {part[w, 1], part[w, 2], part[w, 3]}/Sqrt[part[w, 1]^2 + part[w, 2]^2 + part[w, 3]^2];
      A = {
        {
         part[u, 1] part[v, 1] + part[w, 1]^2 + (part[u, 3] part[w, 2] - part[u, 2] part[w, 3]) (part[v, 3] part[w, 2] - part[v, 2] part[w, 3]),
         part[u, 2] part[v, 1] + part[w, 1] part[w, 2] + (-(part[u, 3] part[w, 1]) + part[u, 1] part[w, 3]) (part[v, 3] part[w, 2] - part[v, 2] part[w, 3]), 
         part[u, 3] part[v, 1] + part[w, 1] part[w, 3] + (part[u, 2] part[w, 1] - part[u, 1] part[w, 2]) (part[v, 3] part[w, 2] - part[v, 2] part[w, 3])
         }, {
         part[u, 1] part[v, 2] + part[w, 1] part[w, 2] + (part[u, 3] part[w, 2] - part[u, 2] part[w, 3]) (-(part[v, 3] part[w, 1]) + part[v, 1] part[w, 3]),
         part[u, 2] part[v, 2] + part[w, 2]^2 + (-(part[u, 3] part[w, 1]) + part[u, 1] part[w, 3]) (-(part[v, 3] part[w, 1]) + part[v, 1] part[w, 3]),
         part[u, 3] part[v, 2] + part[w, 2] part[w, 3] + (part[u, 2] part[w, 1] - part[u, 1] part[w, 2]) (-(part[v, 3] part[w, 1]) + part[v, 1] part[w, 3])
         }, {
         part[u, 1] part[v, 3] + part[w, 1] part[w, 3] + (part[v, 2] part[w, 1] - part[v, 1] part[w, 2]) (part[u, 3] part[w, 2] - part[u, 2] part[w, 3]),
         part[u, 2] part[v, 3] + part[w, 2] part[w, 3] + (part[v, 2] part[w, 1] - part[v, 1] part[w, 2]) (-(part[u, 3] part[w, 1]) + part[u, 1] part[w, 3]),
         part[u, 3] part[v, 3] + (part[u, 2] part[w, 1] - part[u, 1] part[w, 2]) (part[v, 2] part[w, 1] - part[v, 1] part[w, 2]) + part[w, 3]^2
         }
        };
      Join[{part[U, 1] + part[U, 2]}, {part[U, i]}, scale1 U[[3 ;; 5]].A, {scale2 part[U, 6]}], {i, 3, 5}]], 
    CompilationTarget -> "C",
    RuntimeAttributes -> {Listable},
    Parallelization -> True,
    RuntimeOptions -> "Speed"
    ]
   ];
iteration2[data_] := Flatten[citeration2[data], 1];

result2 = NestList[iteration2, data0, 6]; // AbsoluteTiming
Max[Abs[result2 - result]]

(* {0.001042, Null} *)
(* 1.33227*10^-15 *)

 result3 = NestList[iteration2, data0, 9]; // AbsoluteTiming
 Graphics3D[
  Flatten[plot /@ Flatten[result3, 1]], 
  Lighting -> "Neutral", Background -> Black, 
  Boxed -> False, 
  SphericalRegion -> True
 ]

(* {0.018179, Null} *)

The slow part is the rendering by Mathematica, though...

Answer 2

Slow version

Clear["`*"];

s=1./GoldenRatio;
thickness = 0.15;

next[{a_,b_}]:=Table[{a,b}//TranslationTransform[b-a]//
  RotationTransform[Pi/4,{Cos[2k Pi/3],Sin[2k Pi/3],0},b]//
    ScalingTransform[{1,1,1}s,b],{k,3}];

n=5;

pts=NestList[Join@@next/@#&,N@{{{0,0,0},{0,0,1}}},n];//AbsoluteTiming

Graphics3D[{Tube[Join@@pts,0.02]}]

Graphics3D[MapIndexed[With[{id=#2[[1]]},{ColorData["Rainbow",1-id/10],
  MapIndexed[Sphere[#,t=#2[[1]]/10; k=thickness s^id;k(1-t) + t  k s]&,
    Subdivide[#[[1]],#[[2]],9]]}]&,pts,{2}]]

Faster version

Clear["`*"];

s=1./GoldenRatio;
thickness = 0.15;

next=Table[{a,b}//TranslationTransform[b-a]//
 RotationTransform[Pi/4,{Cos[2k Pi/3],Sin[2k Pi/3],0},b]//
   ScalingTransform[{1,1,1}s,b],{k,3}]/.
     Thread[{a,b}->Table[Indexed[A,{x,y}],{x,2},{y,3}]]/.
        expr_:>Compile[{{A,_Real,2}},
           expr,RuntimeAttributes->{Listable} ];

n=9;

pts=MapIndexed[Flatten[#,#2[[1]]-1]&,NestList[next,N@{{{0,0,0},{0,0,1}}},n]];//AbsoluteTiming
(* {0.0085782, Null} *)