On fractal graphics

I was recently browsing, just out of curiosity, an e-book by one of Matlab founder, Cleve Moller called “Experimenting with Matlab“. I recommend it to all new and experienced users of Matlab. This is quite an interesting reading to get started.

On chapter 6, I came across this extremely nice example of fractals in Matlab called Barnsley fern. Please look at the attached picture. This is quite a beautiful example of a fractal that mimic nature quite nicely. Cleve provides a code to make this fractal in Matlab and see it being calculated right under your eyes. Beautifull!

But then, I got frustrated.

This was too slow.

On my machine, it was taking minutes to calculate and display. So I got bored and looked in more details in the code. I launched the profiler and to my surprise, there was something in there that I consider a very bad practice of Matlab, right there, from the very first user of Matlab. I am sure Cleve Moller intention was to show a clear code, not an optimized fast routine. I mean I have eternal respect for one of the writer of the extremely optimized LINPACK routines that are (or were) at the core of Matrix algebra in Matlab. But this is a perfect introduction for me to :

- Emphasize to ALWAYS use the profiler. Even the best of all of us can forget to optimize its code.

- Be carefull when you call graphic display as we will see.

Ok, this is the code, I got from the e-book :

function fern
%FERN  MATLAB implementation of the Fractal Fern.
%   Michael Barnsley, Fractals Everywhere, Academic Press, 1993.
%   This version runs forever, or until stop is toggled.
%   See also FINITEFERN.

shg
clf reset
set(gcf,'color','white','menubar','none', ...
'numbertitle','off','name','Fractal Fern')
x = [.5; .5];
h = plot(x(1),x(2),'.');
darkgreen = [0 2/3 0];
set(h,'markersize',1,'color',darkgreen,'erasemode','none');
axis([-3 3 0 10])
axis off
stop = uicontrol('style','toggle','string','stop', ...
'background','white');
drawnow
p  = [ .85  .92  .99  1.00];
A1 = [ .85  .04; -.04  .85];  b1 = [0; 1.6];
A2 = [ .20 -.26;  .23  .22];  b2 = [0; 1.6];
A3 = [-.15  .28;  .26  .24];  b3 = [0; .44];
A4 = [  0    0 ;   0   .16];
cnt = 1;
tic
while ~get(stop,'value')
r = rand;
if r < p(1)
x = A1*x + b1;
elseif r < p(2)
x = A2*x + b2;
elseif r < p(3)
x = A3*x + b3;
else
x = A4*x;
end
set(h,'xdata',x(1),'ydata',x(2));
drawnow
cnt = cnt + 1;
end

t = toc;
s = sprintf('%8.0f points in %6.3f seconds',cnt,t);
text(-1.5,-0.5,s,'fontweight','bold');
set(stop,'style','pushbutton','string','close','callback','close(gcf)')

The very bad practice in this code is to call the figure window from the while loop at EVERY iteration AND to update ONLY one point at each time. It is even more ridiculous because there are tic and toc calls in the functions as if the timing of the calculations was checked.
All the calls to windows in Matlab are often slow because they require many processes to be done to update the display. So you don’t want to make this call slow down linear algebra calculations which are very very fast (Thanks to Cleve among others).
So the culprits are drawnow and the call to set within the loop. One should vectorize these calls as you would do for any other calculations in matlab.
So I updated the code to this :

function fernfaster

%FERN  MATLAB implementation of the Fractal Fern.
%   Michael Barnsley, Fractals Everywhere, Academic Press, 1993.
%   This version runs forever, or until stop is toggled.
%   See also FINITEFERN.
shg
clf reset
set(gcf,'color','white','menubar','none', ...
'numbertitle','off','name','Fractal Fern')
x = [.5; .5];
h = plot(x(1),x(2),'.');
darkgreen = [0 2/3 0];
set(h,'markersize',1,'color',darkgreen,'erasemode','none');
axis([-3 3 0 10])
axis off
stop = uicontrol('style','toggle','string','stop', ...
'background','white');
drawnow
p  = [ .85  .92  .99  1.00];
A1 = [ .85  .04; -.04  .85];  b1 = [0; 1.6];
A2 = [ .20 -.26;  .23  .22];  b2 = [0; 1.6];
A3 = [-.15  .28;  .26  .24];  b3 = [0; .44];
A4 = [  0    0 ;   0   .16];
cnt = 1;
tic
NumberPlot=2000;
y=zeros(2,NumberPlot);
y(:,1)=x;
k=2;
while ~get(stop,'value')
r = rand;
if r < p(1)
y(:,k) = A1*y(:,k-1) + b1;
elseif r < p(2)
y(:,k) = A2*y(:,k-1) + b2;
elseif r < p(3)
y(:,k) = A3*y(:,k-1) + b3;
else
y(:,k) = A4*y(:,k-1);
end
if k==NumberPlot
CurrentX=get(h,'xdata');
CurrentY=get(h,'ydata');
set(h,'xdata',[CurrentX y(1,2:end)],'ydata',[CurrentY y(2,2:end)]);
k=2;
y(:,1)=y(:,NumberPlot);
drawnow
else
k=k+1;
end
cnt = cnt + 1;
end
t = toc;
s = sprintf('%8.0f points in %6.3f seconds',cnt,t);
text(-1.5,-0.5,s,'fontweight','bold');
set(stop,'style','pushbutton','string','close','callback','close(gcf)')

Note that now, new points are only sent to the window every 2000 points. The code runs much much faster now. I also corrected a small bug. When pressing stop, the window was going back to the first set of points. This was because, the points were not added to the previous xdata but just replacing them. Somehow the display was correct during the processing but not when STOP was pressed.
This is equivalent to problem I had with the waitbar in Matlab. You MUST limit your calls to figures to the minimum necessary, especially within loops.