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.

This entry was posted in Intermediate, Optimizing your code. Bookmark the permalink.

12 Responses to On fractal graphics

  1. This is beautiful!
    Just a very small note: the sign “less than” in the code is shown as the HTML entity, so copy-pasting it gives an error.

  2. [url=http://www.okakaku.com/brand-1-copy-3-cheap-0-max0-attr0.html]確かに気まぐれなエマニュエル・ブーシェの合併症についての1つの時計のように膨大な量があります。今でも、エマニュエル・ブーシェの合併症の2つの作品であるので、私は我々がブーシェ氏の心から期待したいことがたくさんあると思います。[/url]

  3. [url=http://www.bestevance.com/rolex/day-date/index.htm]2012年に戻って、ボールの腕時計は、彼らがビーエムダブリュー車とのコラボレーションで腕時計を製造するという興味深い発表と私はそれはかなりクールだと思った。腕時計のコレクションのためにボールがデビューして、後で行った実際のここで。 ルイヴィトンスーパーコピー 2015年の間、ボールを静かにリリースされた1つの新しい部分をボールにtimetrekkerコレクション青と黒いダイヤルバージョンの両方を含むビーエムダブリュー腕時計と家族のために、鋼の金属ブレスレットやゴムひものどちらかを見て。[/url]

  4. [url=http://www.eevance.com/tokei/ap]1929年アメリカ株式市場の大暴落し、その影響を受けて、自動車製品の販売量を増やすため、アメリカ発起した運動場の技術進歩として、機能が美学さほど重要でない。チュードルスーパーコピーこの運動の後に対してタブ界石英危機に創後捲土重来を生じた啓発を受け、石英危機後のタブ商より重視機械時計はデザイン。この傾向は、90年代末と21世紀初頭の腕時計の作品の中で現れて。[/url]

  5. [url=http://www.gginza.com/%E3%82%A2%E3%83%90%E3%82%A6%E3%83%88/item_1.html]超人気旺店全新登場人気ブランド品★2016年全新登場★揃っている。ロレックス腕時計、オメガ、ブルガリ、カルティエIWC腕時計、シャネルコピーバッグ、財布腕時計、バネライ、ベル&ロス、ゴルム腕時計チュードル腕時計、ルイヴィトン、グッチバッグ財布、コピーミュウミュウバッグ財布、エルメス、ブラダバッグ、財布等弊社は「信用第一」をモットーにお客様にご満足頂けるよう、送料は無料です(日本全国)! ご注文を期待しています!下記の連絡先までお問い合わせください。是非ご覧ください[/url]

  6. [url=http://www.eevance.com/News/6fe97759aa27a0c9.html]ブルガリスーパーコピー 紹介のホリデースペシャルとして、私は1月4日、2016年を通して私の新しいマーク・カーソンkaラスポーツ腕時計は20 %オフで現在提供しています。価格は、色のあなたの選択において2インフィニティ・ストラップを含みます(革のひもは非常にわずかにより多くのコスト)。販売価格1040ドルで始まります。[/url]

  7. [url=http://www.bagkakaku.com/vuitton_bag/2/N41541.html]紹介によると、赫柏林このブランドはいま家族の第3世代の人が経営、赫柏林表も外観デザインや品質の上で、すべて持ってフランスロマンチックと繊細な性格。設計上のオリジナリティも比較的地味、線は優雅な。このブランドの品質について、肖晓もまた、「彼らの品質も守って家族の使命感、利益と家族の栄誉感の選択の中で彼らに後者より。」[/url]

  8. [url=http://www.newkakaku.com/cb11.htm]ボル表(アジア)上海で冠亚時計城の南京西路店で開かれている2009バーゼル鑑賞サロンでは、この中国の消費者に正式に推薦その年度腕時計の新しいデザイン。エルメススーパーコピー今回の露出計5つのシリーズは、2009年にリリースした新しい係デザインだけでなく、全シリーズの汽灯技術普及を続け、更に発揮2009年国際腕時計の設計風潮。[/url]

  9. [url=http://www.eevance.com/News/9d4c2f636f067f89.html]はHamiltonとリーヴァをコラボした「リーヴァtimed by Hamilton」の新シリーズ表項、双方の密接な協力と厳しく要求で、ついには2007年後半に登場。「リーヴァtimed by Hamilton」が最初からのスケッチデザイン1本の表の完成を作って、全行程表にスイス専門技術を懐かしい雰囲気と現代モード感を完璧に解け合って、アイデアの元素を最大限に発揮し。「リーヴァtimed by Hamilton」がそれぞれ出し男ものの時計とパンパン項、一項ごとに表使うのはトップと豪勢な材質をHamilton契合リーヴァや製品の品質に完璧を求めて。スーパーコピー時計同じように無数の嵐を経験して試練のリーヴァヨット世界チャンピオンとしていつも伴っ素手で潜水Pierre Frolla深く潜行のHamilton腕時計、「リーヴァtimed by Hamilton」はリーヴァとHamilton不敵厳しい環境の極緻の傑作。[/url]

  10. [url=http://www.gowatchs.com/brand-214.html]私は最終的にl-evolutionブランパン腕時計の家族はすべて何についての質問に答えますと思いますがそのl-evolutionコレクションの中で真新しい高い合併症は時計をリリースしましたか?」2015年の間、我々はブランパンl-evolutionトゥールビヨンカルーセルその社内で作られた口径2322v2運動の両方の回転木馬とトゥールビヨンは、ダイヤルしましたが、過去にブランパンは2013年により保守的な時計で何かに見えたが含まれます。 ブライトリング スーパーコピー 今、同じ概念と改訂版の動きが存在する現代のスタイルを見る最後の超かもしれないという点でエキゾチックなブランパンl-evolutionウォッチコレクションです。[/url]

  11. [url=http://www.ooowatch.com/tokei/chopard/index.html]およそ1年前、私は私のマーク・カーソンkaラ・スポーツ腕時計はどうなるのだろう、ダイヤルと手を設計したが、特定の思想は当時のストラップに与えられました。私は、バーゼル2015年から帰ったとき、私は座った一つのゴールの設計は、高セキュリティストラップシステムを利用することができたのは、曲がった金属片の「私のマーク・カーソンkaラの場合にnatoストラップのより多くの保安を提供するツールと特徴的なマッチョの観察を提供する。[/url]

  12. [url=http://www.ooobag.com/wallet/louisvuitton/index_11.html]スーパーコピー商品N級当ブランドコピー通販サイトはブランド財布、ブランドバッグのスーパーコピーを専門に 扱っています。コピーブランド 代引き,ルイヴィトンコピー 日本での通販オンラインショップ販売を含め 開業6年の実績を持ち、ご安心してご利用下さい。ブランド N級 ブランドコピー販売 弊社はルイヴィトン、グッチ、ロレックス、オメガ等世界 有名なブランド コピー品を販売しております。スーパーコピーのネット通販最安値が探せるサイトですスーパーコピー専門の販売ショップです。ご安心してお買い物をお楽しみく[/url]

Leave a Reply

Your email address will not be published. Required fields are marked *