my_sierpinski(N,P)

function M = my_sierpinski(N,P)

% 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.
for i=1:max(max(P))*1000

    % 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;
end