function p=randomSurfaceDistribution(v,t,number)
%RANDOMSURFACEDITRIBUTION(vertex,triangle,numberOfPoints) creates vector of (number) points on the
%surface of 3d-shape. the points created are equally distributed over the surface, they are not
%necessarily vertices of the mesh defining the original shape.
%input: array of vertices, array of triangles as given by result of 'read3ds', i.e. the arrays
%       used are the elements of the cell-array-result of read3ds (see example)
%output: (number) x 3 array of points {(x1 y1 z1];...;[xn yn zn]}
%---
%example: plot 1000 points on surface of 3d-object 'file.3ds':
% [vcell,tcell]=read3ds('file.3ds',1);
% p=randomsurfacedistribution(vcell{1},tcell{1},1000);
% plot3(p(:,1),p(:,2),p(:,3),'.');
%---


%idea from 
%'osada, funkhauser et al.: matching 3d models with shape distributions'

%init: create vector of cumulated sorted triangle-areas
s=areaOfTriangles(v,t); %determine area of each triangle
[s,s_index]=sort(s);      %sort
s=cumsum(s);            %cumulate sorted area-values
s=s/s(length(s));       %normalize

%create vector of (number) sorted random numbers
r=rand(number,1);
r=sort(r);

%index triangles by random vector: search max. cumulated area values being
%less or equal than random value
%assuming that rand() produces an equal distribution, the probability of choosing
%a triangle is proportional to its area
index=1;                    %index to s
triangle=zeros(number,1);
for i=1:number
    while r(i)>s(index)
        index=index+1;
    end
    triangle(i)=s_index(index);
end

%now choose a point inside each indexed triangle (one point for each entry)
p=zeros(number,3);
for i=1:number
    tr=t(triangle(i),:);
    a=v(tr(1),:);
    b=v(tr(2),:);
    c=v(tr(3),:);
    r1=rand(1);
    r2=rand(1);
    p(i,:)=(1-sqrt(r1))*a + sqrt(r1)*(1-r2)*b + sqrt(r1)*r2*c;
end

%bye
return


%------------------------------
function s=areaOfTriangles(v,t)

s=zeros(length(t),1);
for i=1:length(t)
    a=v(t(i,1),:)-v(t(i,2),:);
    b=v(t(i,3),:)-v(t(i,2),:);
    s(i)=norm(cross(a,b))/2;
end

return









