function M = my_sierpinski(N,P)

h = figure;
clf reset
xlim([0 1000])
ylim([0 1000])
set(gcf,'color','white','menubar','none', 'numbertitle','off','name','Fractal Sierpinski')

stop = uicontrol('style','toggle','string','stop', 'background','white');
set(stop,'style','pushbutton','string','close','callback','close(gcf)')

% let's create zero matrix N x N
M = zeros(N); 

% Mark each point in matrix M: Change 0 to 1
M(P(1,1), P(1,2)) = 1; 
M(P(2,1), P(2,2)) = 1;
M(P(3,1), P(3,2)) = 1;

% c is current point: reference point
c = P(floor(rand*3)+1 ,:);

% the limitation of i is up to you,
% I thought multiplying max of coordinates by 1000
% would create enough points for the Sierpinski triangle.
% You can reduce it if you like.
hold on;
%i=1;
% i<=max(max(P))*1000

while ~get(stop,'value')

    % for each iteration, create random number between 3 and 1
    % r will help to chose randomly one of the corner points
    r = floor(rand*3)+1; 
    
    % update reference point:
    % Calculate the half distance between randomly chosen corner point
    % and the current point, then new reference point will be this created
    % one.
    c = [floor((c(1) + P(r,1)) / 2), floor((c(2) + P(r,2)) / 2)];
  
    % mark new created point, reference. 0 to 1.
    M(c(1), c(2)) = 1;
    plot(c(1),c(2),'p');
    pause(0.01)
end

hold off;