function RandomMondrian(a,b,L,W,nSubdivides)
% Displays an L-by-W random Mondrian with lower left corner (a,b).
% nSubdivides is a nonegative integer indicating the maximum
% number of subdivisions allowed

pHalt = .2; % Early termination probability. Higher values induce
% bigger tiles.
Bwidth = 2; % Width of tile borders.
f = 2/3; % 0<f<1 The greater the value the more variation in tile size

%plot([a a+L a+L a a],[b b b+W b+W b]);
%plot([a a+L],[b+W/2 b+W/2])
%plot([a+L/2 a+L/2],[b b+W])
% Color map...

MyColors = [ 1.0 0.0 0.0;...
0.0 1.0 0.0;...
0.0 0.0 1.0;...
1.0 0.0 1.0;...
0.0 1.0 1.0;...
1.0 1.0 0.0;...
0.5 0.5 0.0;...
0.5 0.5 0.5;...
0.0 0.5 0.5;...
0.5 0.0 0.5;...
0.8 0.8 0.8];
[n,m] = size(MyColors);

if nSubdivides == 0 || (rand < pHalt && nSubdivides<=2)
% No more subdivisions.
% Choose a random color, fill the rectangle...
k = ceil(rand*n);
fill([a a+L a+L a],[b b b+W b+W],MyColors(k,:),'linewidth',Bwidth)
else
% Pick a random partition point from the middle third of the rectangle...
%c = a + (1+2*rand)*L/4;
%d = b + (1+2*rand)*W/4;
c = a+L*(((1-f)/2)+f*rand());
d = b+W*(((1-f)/2)+f*rand());
% Process each of the four subrectangles...
RandomMondrian(a,b,c-a ,d-b , nSubdivides-1);
RandomMondrian(a,d,c-a ,b+W-d, nSubdivides-1);
RandomMondrian(c,b,a+L-c,d-b , nSubdivides-1);
RandomMondrian(c,d,a+L-c,b+W-d, nSubdivides-1);
end
shg
pause(.05)