N=10000;
p=3;
A=zeros(2,2,3);
c=0.255;
r=0.75;
q=0.625;
phi=-%pi/8;
psi=%pi/5;

A(:,:,1)=[0 0 ; 0 c];
A(:,:,2)=[r*cos(phi) -r*sin(phi) ; r*sin(phi) r*cos(phi)];
A(:,:,3)=[q*cos(psi) -r*sin(psi) ; q*sin(psi) r*cos(psi)];

b=zeros(2,1,3);
b(:,:,1)=[1/2 0];
b(:,:,2)=[1/2-r/2*cos(phi) c-r/2*sin(phi)];
b(:,:,3)=[1/2-q/2*cos(psi) 3*c/5-q/2*sin(psi)];

x=zeros(2,1);
for i=1:N
    u=floor(1+p*rand(1,1,'def'));
    x=A(:,:,u)*x+b(:,:,u);
    plot2d(x(1),x(2),0,rect=[0,0,1,1])
end