> restart:

> h:=proc(c,r,m)
RETURN(c+r*(m-c))
end:

> s:=proc(c,r,t,m) local xc,yc,xm,ym:
xc:=op(1,c):yc:=op(2,c):
xm:=op(1,m):ym:=op(2,m):
c+[r*(cos(t)*(xm-xc)-sin(t)*(ym-yc)),r*(sin(t)*(xm-xc)+cos(t)*(ym-yc))]:
end:

> transfo:=proc(l,theta)
local A,B,r1,r2:
A:=op(1,l):B:=op(2,l):r1:=1/4/cos(theta)**2:r2:=1/2/cos(theta):
[A,h(A,r1,B),s(A,r2,theta,B),h(B,r1,A),B]:
end:

> m:=8:
A:=[0,0]:B:=[1,0]:
Digits:=3:for k from 0 to m do
l:=[A,B]:theta:=evalf(Pi/6+k*Pi/12/m):
n:=5:
for i to n do
ll:=[A]:
for j to 4^(i-1) do
lp:=[op(j,l),op(j+1,l)]:
ll:=[op(ll),op(subsop(1=NULL,transfo(lp,theta)))]:
od:
l:=ll:
od:
q:=floor(nops(l)/2)+1:
l.k:=[seq(op(map(t->s(A,1,-(Pi-4*theta),[t[1],-t[2]]),l[1..q]))[q-i+1],i=1..q),op(l),seq(op(map(t->s(B,1,Pi-4*theta,[t[1],-t[2]]),l[q..nops(l)]))[q-i+1],i=1..q)]
od:

> display([seq(plot(l.j,thickness=2),j=0..m),seq(plot(l.(m-j),thickness=2),j=1..m-1)],insequence=true,axes=none,scaling=constrained);