# How to draw a Pythagoras tree like this

I want to draw a Pythagoras tree like the one below:

But I can only do this now, and still can't do random colors.

\[ScriptCapitalD] = DiagonalMatrix[{1., 1.}];
\[ScriptCapitalD]1 = DiagonalMatrix[{1., 1.}/(Sqrt[2])];
ifs = {Composition[AffineTransform[{\[ScriptCapitalD], {0, 0}}],
RotationTransform[0]],
Composition[AffineTransform[{\[ScriptCapitalD]1, {0, 1.}}],
RotationTransform[Pi/4]],
Composition[AffineTransform[{\[ScriptCapitalD]1, {1./2, 3./2}} ],
RotationTransform[-Pi/4]]};
g = Graphics[Rectangle[]];
Transform[0] := g;
Transform[n_Integer?Positive] := Transform[n] =
Graphics @
Table[GeometricTransformation[Transform[n - 1] // First,
ifs[[i]]], {i, Length[ifs]}];

Table[Transform[n], {n, 0, 8}]


Here is one possibility:

branch[Polygon[{u_, v_, r_, s_}], q_?NumericQ] := Module[{a = Norm[r - s], b, h, v1, v2},
b = a Sqrt[1 - q^2];
h = r q^2 + s (1 - q^2) + b q Normalize[Cross[r - s]];
v1 = a q Cross[Normalize[h - s]]; v2 = -b Cross[Normalize[h - r]];
{Polygon[{s, h, h + v1, s + v1}], Polygon[{h, r, r + v2, h + v2}]}]

Block[{a = 1., q = 0.5, n = 9},
Graphics[MapIndexed[{ColorData[63] @@ #2, #1} &,
NestList[With[{qq = (q = Sqrt[1 - q^2])},
# /. p_Polygon :> branch[p, qq]] &,
Polygon[{{0, 0}, {a, 0}, {a, a}, {0, a}}], n]],
Background -> Black]]


Clear["*"];
next = Table[Indexed[A,{x,y}],{x,4},{y,2}]/.
v:{a_,b_,c_,d_}:>
Compile[{{A,_Real,2},t},
{
ScalingTransform[{1,1}Cos[t],d]@RotationTransform[t,d]@TranslationTransform[d-a]@v,
ScalingTransform[{1,1}Sin[t],c]@RotationTransform[t-Pi/2,c]@TranslationTransform[c-b]@v
}//Evaluate,
RuntimeAttributes->{Listable}
];

n = 15;
t = Pi/3;
pts = NestList[next[#,t=Pi/2-t]&, N@{{{0,0},{1,0},{1,1},{0,1}}},n];
poly = MapIndexed[With[{i=Tr@#2},{ColorData["Rainbow",i/(n+1)],Polygon@Flatten[#,i-1]}]&, pts];
Graphics[{Antialiasing->True, poly}, Background->Black,ImageSize->600]


generalization to 3D is also easy

Clear["*"];
next=Table[Indexed[A,{x,y}],{x,8},{y,3}]/.
v:{a1_,b1_,c1_,d1_,a2_,b2_,c2_,d2_}:>
Compile[{{A,_Real,2},t},
{

ScalingTransform[{1,1,1}Cos[t],(a2+d2)/2]@RotationTransform[t,{0,-1,0},(a2+d2)/2]@TranslationTransform[a2-a1]@v,
ScalingTransform[{1,1,1}Sin[t],(b2+c2)/2]@RotationTransform[Pi/2-t,{0,1,0},(b2+c2)/2]@TranslationTransform[a2-a1]@v
}//Evaluate,
RuntimeAttributes->{Listable}
];

n=10;
t=Pi/3;
pts=NestList[next[#,t=Pi/2-t]&,N@{{{0,0,0},{1,0,0},{1,1,0},{0,1,0},{0,0,1},{1,0,1},{1,1,1},{0,1,1}}},n];
poly=MapIndexed[With[{i=Tr@#2},{Blend[{Cyan,Green,Yellow},i/(n+1)],Hexahedron@Flatten[#,i-1]}]&,pts];
Graphics3D[{EdgeForm[],poly},Boxed->False,Background->RGBColor[{0,137,188}/255],ImageSize->600,Lighting->"Neutral"]


• How do you efficiently create step-by-step animation? Jan 31, 2020 at 0:13