MATLAB Compiler

A Simple Example - The Sierpinski Gasket

Consider an M-file function called gasket.m.

function theImage = gasket(numPoints)
%GASKET An image of a Sierpinski Gasket.
%   IM = GASKET(NUMPOINTS)
%
%   Example:
%   x = gasket(50000);
%   imagesc(x);colormap([1 1 1;0 0 0]);
%   axis equal tight

%   Copyright (c) 1984-98 by The MathWorks, Inc
%   $Revision: 1.1 $  $Date: 2000/09/08 18:42:05 $

theImage = zeros(1000,1000);

corners = [866 1;1 500;866 1000];
startPoint = [866 1];
theRand = rand(numPoints,1);
theRand = ceil(theRand*3);

for i=1:numPoints
   startPoint = floor((corners(theRand(i),:)+startPoint)/2);
   theImage(startPoint(1),startPoint(2)) = 1;
end

How the Function Works

This function determines the coordinates of a Sierpinski Gasket using an Iterated Function System algorithm. The function starts with three points that define a triangle, and starting at one of these points, chooses one of the remaining points at random. A dot is placed at the midpoint of these two points. From the new point, a dot is placed at the midpoint between the new point and a point randomly selected from the original points. This process continues and eventually leads to an approximation of a curve.

The curve can be graphed in many ways. Sierpinski's method is:

• Start with a triangle and from it remove a triangle that is one-half the height of the original and inverted. This leaves three triangles.
• From each of the remaining three triangles, remove a triangle that is one-fourth the height of these new triangles and inverted. This leaves nine triangles.
• The process continues and at infinity the surface area becomes zero and the length of the curve is infinite.

gasket.m is a good candidate for compilation because it contains a loop. The overhead of the for loop command is relatively high compared to the cost of the loop body. M-file programmers usually try to avoid loops containing scalar operations because loops run relatively slowly under the MATLAB interpreter.

To achieve a reasonable approximation of the Sierpinski Gasket, set the number of points to 50,000. To compute the coordinates and time the computation, you can use

tic; x = gasket(50000); toc

To display the figure, you can use

imagesc(x); colormap([1 1 1;0 0 0]);
axis equal tight

Getting Started with MEX-Files

Invoking the M-File