Creating Sierpinski gasket fractal with MATLAB

Here is a simple code to create Sierpinski fractals, you can refer to "Fractal Image Compression Using Iterated Transforms: Application to DTED" for more detail

here is a function to apply single iteration transform of W=w1 U w2 U w3, where U represents set union.

function [ tR ] = Wn(S)
%Creating Sierpinski gasket fractal

[r c]=size(S);
blank=ones(r,c);
rsz=[0.5 0;
     0 0.5];
d1=[0; 
    0];
d2=[0;
    fix(r/2)];
d3=[fix(r/2);
    0];
x=repmat(1:r,[r 1]);y=x';x=x(:);
y=y(:);
xy=[x y]';

d1=repmat(d1,[1 size(xy,2)]);
newXY1=rsz*xy+d1;
newXY1=fix(newXY1);newXY1(newXY1==0)=1;
blank(newXY1(1,:),newXY1(2,:))=S(xy(1,:),xy(2,:));

d2=repmat(d2,[1 size(xy,2)]);
newXY2=rsz*xy+d2;
newXY2=fix(newXY2);newXY2(newXY2==0)=1;
blank(newXY2(1,:),newXY2(2,:))=S(xy(1,:),xy(2,:));

d3=repmat(d3,[1 size(xy,2)]);
newXY3=rsz*xy+d3;
newXY3=fix(newXY3);newXY3(newXY3==0)=1;
blank(newXY3(1,:),newXY3(2,:))=S(xy(1,:),xy(2,:));

tR=blank;
end

The above function can be called to apply only the first iteration, you can call the function as many times as you wish using the following code

R=zeros(50,50);
for i=1:4
    R=Wn(R);
end

figure, imshow(R,[]);

which creates a fractal of size 50 x 50, but with only 4 iterations, the output is shown in the following figure