Square fractal generation in tikz

I'm trying to recreate this exact diagram using tikz. So far I have been attempting to use nested decorations in tikz since I used these for other fractal constructions. However, these constructions were all either decorations pre-defined by tikz like the Koch curve or ones I found solutions to on stack exchange, like the Sierpinski triangle.

I've been looking a lot into defining my own decorations but it seems like quite a complicated process for a tikz newbie, and haven't found any examples too similar to what I'm trying to do. I know it will also be possible using lindemayer systems, but only understand how to use them for line constructions.

If it's any help, to my mind it seems the simplest way to do this would be to set the square as the initial shape with origin at the bottom left, then do scale by 1/4 for the bottom left square, scale by 1/4 then translate upwards for the top left square etc. and then set the new shape to replace the initial, ready for the next iteration.

Any help would be much appreciated!

Here's a way with a Lindenmayer system. For orders above 5, compile with LuaLaTeX.

% \RequirePackage{luatex85} % Only for LuaLaTeX and standalone class
\documentclass[varwidth,border=5]{standalone}
\usepackage{tikz}
\usetikzlibrary{lindenmayersystems}
\pgfdeclarelindenmayersystem{square fractal}{%
\symbol{S}{\pgflsystemstep=0.5\pgflsystemstep}
\symbol{A}{\pgftransformshift%
{\pgfqpoint{0.75\pgflsystemstep}{0.75\pgflsystemstep}}}
\symbol{R}{\pgftransformrotate{90}}
\symbol{Q}{%
\pgfpathrectangle{\pgfqpoint{-0.5\pgflsystemstep}{-0.5\pgflsystemstep}}%
{\pgfqpoint{\pgflsystemstep}{\pgflsystemstep}}%
}
\rule{Q -> [SQ[ASQ][RASQ][RRASQ][RRRASQ]]}
}
\begin{document}
\foreach\i in {0,...,5}{%
\tikz\fill [l-system={square fractal, step=5cm, axiom=Q, order=\i}]
lindenmayer system;
\ifodd\i\par\bigskip\leavevmode\fi
}
\end{document}


And here's a way with decorations:

\documentclass[varwidth,border=5]{standalone}
\usepackage{tikz}
\usetikzlibrary{decorations}
\pgfdeclaredecoration{square fractal}{start}{
\state{start}[width=0pt,next state=draw]{
\pgfpathmoveto{\pgfpointdecoratedinputsegmentfirst}
}
\state{draw}[width=\pgfdecoratedinputsegmentlength]{
\pgfpointdiff{\pgfpointdecoratedinputsegmentfirst}%
{\pgfpointdecoratedinputsegmentlast}
\pgfgetlastxy\tmpx\tmpy
\pgfmathveclen\tmpx\tmpy
\pgfmathparse{\pgfmathresult/4}%
\let\tmp=\pgfmathresult
\pgfpathlineto{\pgfpoint{\tmp}{+0pt}}
\pgfpathlineto{\pgfpoint{\tmp}{-\tmp}}
\pgfpathlineto{\pgfpoint{3*\tmp}{-\tmp}}
\pgfpathlineto{\pgfpoint{3*\tmp}{+0pt}}
\pgfpathlineto{\pgfpointdecoratedinputsegmentlast}
}
\state{final}{
\pgfpathclose
}
}
\begin{document}
\tikz[decoration=square fractal]
\fill (0,0) rectangle (4,4);
\tikz[decoration=square fractal]
\fill decorate { (0,0) rectangle (4,4) };
\\
\tikz[decoration=square fractal]
\fill decorate { decorate { (0,0) rectangle (4,4) } };
\tikz[decoration=square fractal]
\fill decorate { decorate { decorate { (0,0) rectangle (4,4) } } };
\end{document}


TikZ solution

The black squares of the fractal are generated via an expandable recursion.

\documentclass[tikz]{standalone}

\usepackage{etoolbox}
\makeatletter
\patchcmd{\tikz@@command@path}{=100}{=10000}{}{\errmessage{Patching failed.}}
\makeatother

\makeatletter
\newcommand*{\@SquareFractal}[4]{%
% #1: order
% #2: edge length
% #3: x position of lower left corner
% #4: y position of lower left corner
\ifnum#1=0
(#3,#4)rectangle(\the\dimexpr(#3)+(#2)\relax,\the\dimexpr(#4)+(#2)\relax)%
\expandafter\@gobble
\else
\expandafter\@firstofone
\fi
{
% Middle
\expandafter\@SquareFractal
\expandafter{\the\numexpr(#1)-1\expandafter}%
\expandafter{\the\dimexpr(#2)/2\expandafter}%
\expandafter{\the\dimexpr(#3)+(#2)/4\expandafter}%
\expandafter{\the\dimexpr(#4)+(#2)/4}%
% Bottom left
\expandafter\@SquareFractal
\expandafter{\the\numexpr(#1)-1\expandafter}%
\expandafter{\the\dimexpr(#2)/4}%
{#3}%
{#4}%
% Bottom right
\expandafter\@SquareFractal
\expandafter{\the\numexpr(#1)-1\expandafter}%
\expandafter{\the\dimexpr(#2)/4\expandafter}%
\expandafter{\the\dimexpr(#3)+(#2)*3/4}%
{#4}%
% Top left
\expandafter\@SquareFractal
\expandafter{\the\numexpr(#1)-1\expandafter}%
\expandafter{\the\dimexpr(#2)/4\expandafter}%
\expandafter{\the\dimexpr(#3)\expandafter}%
\expandafter{\the\dimexpr(#4)+(#2)*3/4}%
% Top right
\expandafter\@SquareFractal
\expandafter{\the\numexpr(#1)-1\expandafter}%
\expandafter{\the\dimexpr(#2)/4\expandafter}%
\expandafter{\the\dimexpr(#3)+(#2)*3/4\expandafter}%
\expandafter{\the\dimexpr(#4)+(#2)*3/4}%
}%
}

\newcommand*{\SquareFractal}[2]{%
% #1: order
% #2: edge length
\begingroup
\edef\x{\@SquareFractal{#1}{#2}{0pt}{0pt}}%
\expandafter\tikz\expandafter\fill\x;%
\endgroup
}
\makeatother

\begin{document}
\foreach\i in {0, ..., 5} {\SquareFractal{\i}{\linewidth}}
\end{document}


Since the whole drawing commands are hold in memory, the memory is the limiting factor.

Result for order 5:

IniTeX solution

The following example uses simple rules in iniTeX for drawing the squares to get higher orders without running out of memory.

The maximum dimension in TeX is 16383.99998 pt (\maxdimen). This is (230 - 1) sp (1 pt = 216 sp = 65536 sp). The smallest squares of the next level uses a square edge length of a quarter. Then, it follows that with the smallest square edge length of 1 sp, the largest order is 14, the edge length of the result is then 228 sp.

The example uses either pdfTeX or luaTeX in iniTeX mode (pdftex -ini -etex or luatex -ini). LuaTeX is faster and has less memory restrictions. For comparison, order 8 takes around 45 s with pdfTeX, but 8 s with LuaTeX. Higher orders with LuaTeX:

• Order 10: time is 3 3/4 min, file size is 47 MiB.

• Order 11: time is 33 min, file size is 173 MiB.

At order 12, the computer gave up and I had to reboot.

Example:

\catcode\{=1
\catcode\}=2
\catcode\#=6

\ifx\directlua\undefined
\pdfoutput=1
\pdfminorversion=4
\pdfhorigin=0pt
\pdfvorigin=0pt
\pdfcompresslevel=9
\else
\directlua{%
tex.enableprimitives('', {'outputmode', 'dimexpr', 'numexpr'})
tex.enableprimitives('pdf', {'pagewidth', 'pageheight'})
}
\outputmode=1
\directlua{
pdf.setorigin()
pdf.setminorversion(4)
pdf.setcompresslevel(9)
}
\fi

\dimendef\pagewidth=0
\dimendef\xpos=2

\def\SquareFractal#1#2{%
% #1: order
% #2: minimum edge length
\pagewidth=\dimexpr#2\MulFour#1!\relax
\immediate\write16{* Calculating square fractal of order #1 ...}%
\pdfpagewidth=\pagewidth %
\pdfpageheight=\pagewidth %
\shipout\hbox{%
\xpos=0pt\relax
\SquareFractalRecursiv#1!\pagewidth!0pt!0pt!%
\kern\dimexpr\pagewidth-\xpos\relax
}%
}

\def\MulFour#1!{%
\ifnum#1=0
\else
*4%
\expandafter\MulFour
\the\numexpr#1-1\expandafter!%
\fi
}

\def\SquareFractalRecursiv#1!#2!#3!#4!{%
% #1: order
% #2: edge length
% #3: x position of lower left corner
% #4: y position of lower left corner
\ifnum#1=0 %
\iffalse
\raise#4\hbox to 0pt{%
\kern#3\relax
\vrule width#2height#2\relax
\hss
}%
\else
\ifdim#3=\xpos
\else
\kern\dimexpr#3-\xpos\relax
\fi
\vrule width#2 depth-#4 height\dimexpr#4+#2\relax
\xpos=\dimexpr#3+#2\relax
\fi
\else
% Lower left square
\expandafter\SquareFractalRecursiv
\the\numexpr#1-1\expandafter!%
\the\dimexpr#2/4\expandafter!%
#3!%
#4!%
% Middle square
\expandafter\SquareFractalRecursiv
\the\numexpr#1-1\expandafter!%
\the\dimexpr#2/2\expandafter!%
\the\dimexpr#3+#2/4\expandafter!%
\the\dimexpr#4+#2/4!%
% Lower right square
\expandafter\SquareFractalRecursiv
\the\numexpr#1-1\expandafter!%
\the\dimexpr#2/4\expandafter!%
\the\dimexpr#3+#2*3/4!%
#4!%
% Upper left square
\expandafter\SquareFractalRecursiv
\the\numexpr#1-1\expandafter!%
\the\dimexpr#2/4\expandafter!%
\the\dimexpr#3\expandafter!%
\the\dimexpr#4+#2*3/4!%
% Upper right square
\expandafter\SquareFractalRecursiv
\the\numexpr#1-1\expandafter!%
\the\dimexpr#2/4\expandafter!%
\the\dimexpr#3+#2*3/4\expandafter!%
\the\dimexpr#4+#2*3/4\expandafter!%
\fi
}

% BTW, unit bp instead of pt decreases the output file size
% a bit because of less fractional digits.

% \SquareFractal{<order>}{<length of smallest square>}
% The values of the follwing calls are used in such a way
% that the generated fractals with different orders have
% the same widths and heights.

\SquareFractal{0}{4096pt}
\SquareFractal{1}{1024pt}
\SquareFractal{2}{256pt}
\SquareFractal{3}{64pt}
\SquareFractal{4}{16pt}
\SquareFractal{5}{4pt}
\SquareFractal{6}{1pt}% 65536 sp
\SquareFractal{7}{16384sp}
\SquareFractal{8}{4096sp}
\SquareFractal{9}{1024sp}
\SquareFractal{10}{256sp}
\SquareFractal{11}{64sp}
% \SquareFractal{12}{16sp}
% \SquareFractal{13}{4sp}
% \SquareFractal{14}{1sp}
\end


Result for order 11 (better resolutions are rejected by imgur):

Because of the sheer number of squares, viewing a PDF with higher orders slows down the PDF viewer.

It is therefore more efficient to generate a monochrome bitmap image, e.g. with the smallest squares as squares of 1 x 1 pixel. The image width and height for order 11 is then 222 pixel = 4194304 pixel.

Here is an attempt with MetaPost, for whom it may interest. The recursive macro (square_fractal) at the basis of this program is heavily inspired by this answer to a closely related subject.

vardef square_fractal(expr A, B, n) =
save P; pair P[]; P0 = A; P1 = B;
for i = 1 upto 2:
P[i+1] = P[i-1] rotatedaround (P[i], -90);
endfor;
if n = 0: fill P0 for i = 1 upto 3: -- P[i] endfor -- cycle;
else:
save Q; pair Q[];
for i = 0, 2:
Q[i] = 1/4[P[i],P[i+1]]; Q[i+1] = 3/4[P[i],P[i+1]];
square_fractal(P[i], Q[i], n-1);
square_fractal(Q[i+1], P[i+1], n-1);
endfor;
square_fractal(P0 rotatedaround (Q0, -90), P1 rotatedaround (Q1, 90), n-1); fi
enddef;

beginfig(1);
for n = 0 upto 4:
draw image(square_fractal(origin, (4cm, 0), n)) shifted (n*4.5cm, 0);
endfor;
endfig;

end.


Starting from order 0 (the full square), MetaPost manages an output up to order 6 on my machine. Interestingly enough, order 7 is reached if the previous code is included in a LuaLaTeX program. I don't know the reason why.

Edit Still within LuaLaTeX, and after using floating point numerics (\mplibnumbersystem{double} added right after \usepackage{luamplib}) instead of the defaults fixed point numerics , MetaPost manages to produce the figure at order 9 after 20 minutes. But it nearly freezes my very old laptop (a MacBook Pro from 2008), so I don't dare to go further on it. Maybe I will try it again on a more recent and more powerful computer.

\RequirePackage{luatex85}
\documentclass[border=5mm]{standalone}
\usepackage{luamplib}
\mplibnumbersystem{double}
\begin{document}
\begin{mplibcode}

vardef square_fractal(expr A, B, n) =
save P; pair P[]; P0 = A; P1 = B;
for i = 1 upto 2:
P[i+1] = P[i-1] rotatedaround (P[i], -90);
endfor;
if n = 0: fill P0 for i = 1 upto 3: -- P[i] endfor -- cycle;
else:
save Q; pair Q[];
for i = 0, 2:
Q[i] = 1/4[P[i],P[i+1]]; Q[i+1] = 3/4[P[i],P[i+1]];
square_fractal(P[i], Q[i], n-1);
square_fractal(Q[i+1], P[i+1], n-1);
endfor;
square_fractal(P0 rotatedaround (Q0, -90), P1 rotatedaround (Q1, 90), n-1); fi
enddef;

beginfig(1);
square_fractal(origin, (12cm, 0), 9);
endfig;

\end{mplibcode}
\end{document}


The figure below is the one of order 8. I couldn't manage to produce a PNG version of order 9 because of the near-freezing of my laptop.

Here is a brute force code, far less elegant than Mark Wilbrows code, but perhaps suitable to learn TikZ basics. (There is one slight difference, though, namely in this version the branches extend in all directions except to the one where they came from. This requires the \ifnum statement involving \Veto, which prevents the branches from growing "backwards".)

\documentclass{article}
\usepackage{tikz}
\usetikzlibrary{calc}
\tikzset{blackbox/.style ={rectangle,minimum width=1cm, minimum
height=1cm,scale=#1,draw,fill=black}}

\newcommand{\Fractal}[4][]{
\pgfmathsetmacro{\BoxWidth}{#2 cm}
\node[#1,blackbox=#2] (#3) at #4 {};
\foreach \i in {1,...,4}
{
\pgfmathsetmacro{\HalfWidth}{\BoxWidth/2}
\node[#1,blackbox={#2/2}] (#3-\i) at ($#4+({\i*90-45}:\BoxWidth pt)$) {};
}
}
\begin{document}
\begin{tikzpicture}
\Fractal{1}{0}{(0,0)}
\foreach \j in {1,...,4}
{
\Fractal{0.5}{\j}{({\j*90-45}:2cm)}
\foreach \k in {1,...,4}
{
\pgfmathtruncatemacro{\Veto}{mod(\j-\k+2,4)}
\ifnum\Veto=0\relax
\else
\Fractal{0.25}{\k}{($({\j*90-45}:2cm)+({\k*90-45}:1cm)$)}
\fi
}
}
\end{tikzpicture}
\end{document}


Another alternative with Tikz and recursion.

\documentclass{article}
\usepackage{tikz}
\usetikzlibrary{calc}
\newcommand\DrawFracSquare[4]{{% {Current number}{Side Length}{X}{Y}
\ifnum#1=0
\fill[black] ($(#3,#4)-(#2/2,#2/2)$) rectangle +(#2,#2);
\else
\pgfmathsetmacro\NewNumber{int(#1-1)}
\pgfmathsetmacro\NewSideLength{#2/2}
\edef\NewRec{\noexpand\DrawFracSquare{\NewNumber}{\NewSideLength}{#3}{#4}}
\NewRec
\pgfmathsetmacro\NewSideLength{#2/4}
\pgfmathsetmacro\NewX{#3+3*#2/8}
\pgfmathsetmacro\NewY{#4+3*#2/8}
\edef\NewRec{\noexpand\DrawFracSquare{\NewNumber}{\NewSideLength}{\NewX}{\NewY}}
\NewRec
\pgfmathsetmacro\NewX{#3-3*#2/8}
\pgfmathsetmacro\NewY{#4+3*#2/8}
\edef\NewRec{\noexpand\DrawFracSquare{\NewNumber}{\NewSideLength}{\NewX}{\NewY}}
\NewRec
\pgfmathsetmacro\NewX{#3-3*#2/8}
\pgfmathsetmacro\NewY{#4-3*#2/8}
\edef\NewRec{\noexpand\DrawFracSquare{\NewNumber}{\NewSideLength}{\NewX}{\NewY}}
\NewRec
\pgfmathsetmacro\NewX{#3+3*#2/8}
\pgfmathsetmacro\NewY{#4-3*#2/8}
\edef\NewRec{\noexpand\DrawFracSquare{\NewNumber}{\NewSideLength}{\NewX}{\NewY}}
\NewRec
\fi
}}
\begin{document}
\begin{tikzpicture}
\DrawFracSquare{0}{3}{0}{4}
\DrawFracSquare{1}{3}{4}{4}
\DrawFracSquare{2}{3}{8}{4}
\DrawFracSquare{3}{3}{0}{0}
\DrawFracSquare{4}{3}{4}{0}
\DrawFracSquare{5}{3}{8}{0}
\end{tikzpicture}
\end{document}
`