How to draw a Pythagoras tree like this

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]]

Related link:
https://mp.weixin.qq.com/s/oz1c3crqgC5WMbUFGxLDMQ https://codegolf.stackexchange.com/questions/25043/i-like-pythagorean-trees/210327#210327

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"]