function sierpinski_random(Nstep)
	%% Sierpinski fractal - random points method

	if nargin < 1, Nstep=1500000; end

	% triangle base
	tBase=800;
	
	points=zeros(Nstep,2);

	% triangle vertexes coordinates
	% upside triangle 
	v(1, :)=[1, 1];
	v(2, :)=[tBase, 1];
	v(3, :)=[ ceil(tBase/2), tBase];

	
	% initial point: any point is good as long as it inside triangle
	p= v(1,:);

	for i=1:Nstep
		%place point to the image
		points(i,:) = [ p(1),  p(2) ];

		%choose a random vertex
		n=ceil ( rand(1)*3 );

		% new point
		p= ( ( p + v(n,:) )/2 ); % midpoint between  old one and the vertex
	end

	% % fill the image matrix
	%img=zeros(tBase, tBase);
	%for i=1:Nstep
		% % recall that (y,x) = (1,1) is left upper conner
		%img( ceil(tBase-points(i,2)) , ceil(points(i,1) ) ) =1;
	%end
	%figure()
	%imagesc(img);

	figure()
	plot( points(:,1), points(:,2),   '.');

end

	




