Efficient code for draw a sierpinski triangle with matlab

Question

Following is my code:

function sierpinski(A, B, C, n)
    if n == 0
        patch([A(1), B(1), C(1)], [A(2), B(2), C(2)], [0.0 0.0 0.0]);
    else
       sierpinski(A, (A + B)/2, (A + C)/2, n-1);
       sierpinski(B, (B + A)/2, (B + C)/2, n-1);
       sierpinski(C, (C + A)/2, (C + B)/2, n-1);
end

% sierpinski([0 0], [1 0], [.5 .8], 8)

It's not very effectly. I want to first generating all data then patched, but I don't know how to correctly used. Also, can my code be written use for loops?

Answer 1

Your idea to write one function to generate the data and another to plot it is a good one - it's often a good idea to separate data generation from processing, and processing from output. I would do it something like this:

function out = sierpinski(a, b, c, n)

    if n == 0
        out.xvals = [a(1), b(1), c(1)];
        out.yvals = [a(2), b(2), c(2)];
    else
        out1 = sierpinski(a, (a+b)/2, (a+c)/2, n-1);
        out2 = sierpinski(b, (a+b)/2, (b+c)/2, n-1);
        out3 = sierpinski(c, (a+c)/2, (b+c)/2, n-1);
        out = [out1, out2, out3];
    end

end

This creates a struct of length 3^n, each entry of which contains the coordinates of one of the small triangles in the sierpinski triangle. Your code to plot it might then look like

>> out = sierpinski([0,0], [1,0], [0.5, sqrt(3)/2], 8);
>> figure(); hold on;
>> for i = 1:length(out)
       patch(out(i).xvals, out(i).yvals, 'k');
   end

That crashes on my machine (it seems that Matlab doesn't handle thousands of patches on the same plot very well) but a similar loop which plots one point at the corner of each small triangle.

>> x = [out.xvals];
>> y = [out.yvals];
>> plot(x, y, '.');

which produces this plot

Answer 2

I don't have any code-examples ready, but:

Instead of drawing each triangle as a single patch-object, you could try drawing all triangles in one large patch. Basically you'd only need to concatenate the x- and y-coordinates for each triangle separated by a NaN - this will prevent the patch from drawing lines connecting individual triangles. E.g. the following line produces two separate triangles:

p = patch( [0 0.5 1 0 NaN 2 2.5 3 2 NaN ], [ 0 1 0 0 NaN 2 3 2 2 NaN], 'k')

Mind that, to have a closed triangle you need 4 points per triangle this way, the last point being identical to the first.

Answer 3

EDIT: For each recursion level, you can 'erase' the central triangles, so you have to patch many fewer triangles. For example, at first level, you have three 'up' triangles and only one 'down' triangle. You can path this, instead of the other three. a more compact routine is:

function sierpinski(rec)
  [x, x0] = deal(cat(3, [1 0]', [-1 0]', [0 sqrt(3)]')); 
  for k = 1 : rec x = x(:,:) + x0 * 2 ^ k / 2; 
end 
patch('Faces', reshape(1 : 3 * 3 ^ k, 3, '')', 'Vertices', x(:,:)') 
end

---

So you have to fill much less triangles if...

function sierpinski(rec)
close all
%Main Triangle
hFig=figure;
units=get(hFig,'units');
set(hFig,'units','normalized','outerposition',[0 0 1 1], 'Color', 'white');
set(hFig,'units',units); clear units
hold on
Vx=[0 0.5 1]; Vy=[0 realsqrt(3)/2 0];
fill(Vx,Vy,'b')
%the number of white triangles = sum(3.^(0:1:rec-1))
whitex=NaN(3,sum(3.^(0:1:rec-1))); whitey=whitex; K=1;
for S=1:rec
    [Vx,Vy]=sierpinskisect;
end
fill(whitex,whitey,'w')

function [outX,outY]=sierpinskisect
    %the number of blue triangles = 3^S
    L=size(Vx,1);
    outX=NaN(3*L,3); outY=outX; J=1;
    for I=1:L
        %left blue triangle
        outX(J,:)=[Vx(I,1) mean(Vx(I,(1:2))) mean(Vx(I,([1 3])))];
        outY(J,:)=[Vy(I,1) mean(Vy(I,(1:2))) mean(Vy(I,([1 3])))];
        J=J+1;
        %right blue triangle
        outX(J,:)=[mean(Vx(I,([1 3]))) mean(Vx(I,(2:3))) Vx(I,3)];
        outY(J,:)=[mean(Vy(I,([1 3]))) mean(Vy(I,(2:3))) Vy(I,3)];
        J=J+1;
        %upper blue triangle
        outX(J,:)=[mean(Vx(I,(1:2))) Vx(I,2) mean(Vx(I,(2:3)))];
        outY(J,:)=[mean(Vy(I,(1:2))) Vy(I,2) mean(Vy(I,(2:3)))];
        J=J+1;
        %white triangle
        whitex(:,K)=[outX(J-3,2);outX(J-3,3);outX(J-2,2)];
        whitey(:,K)=[outY(J-3,2);outY(J-3,3);outY(J-2,2)];
        K=K+1;
    end
end

end