# 3D tree in Mathematica?

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

• This and this. – corey979 Aug 21 '17 at 23:54
• @ corey979 Thank you very much for the information, read it with great attention to understand the methods for the generation of this type of graphics. – bullitohappy Aug 22 '17 at 4:49
• In case you don't know, it's common courtesy to link to the corresponding Community question if you ask something here and at Community. (Likewise, if you post there, link to your question here.) – J. M.'s torpor Aug 22 '17 at 10:37
• I want to apologize to everyone for not putting the link that leads to the community as you suggest, I really did not know, but I will not forget, greetings. – bullitohappy Aug 27 '17 at 0:39
• You might be interested in so-called "algorithmic botany" or "computational botany". – Silvia Aug 6 '19 at 8:16

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 = CompileGetElement},
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...

• What is getp? – J. M.'s torpor Aug 22 '17 at 4:03
• Ah. Wait a sec. I am still building on this... – Henrik Schumacher Aug 22 '17 at 4:04
• @J.M. Hm. I guess, it's done. Have a look. – Henrik Schumacher Aug 22 '17 at 4:31
• @HenrikSchumacher Astonishing, the tree is the same as the one on the website, mentioned separately are the improvements you propose. I am learning more and more with these kinds of answers. Thanks for your answer – bullitohappy Aug 22 '17 at 4:57
• You're welcome! Actually, it was quite some fun! – Henrik Schumacher Aug 22 '17 at 5:54

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}]/.
`