function [cd, centroid, samples]=centroidDistance(P,n,show)
%CENTROIDDISTANCE(P,n,<show> centroid distance function of polygon P sampled at n points
%equally spaced along the boundary
%input: Polygon P [x0 y0; x1 y1;..;xn yn]
%       number of sample points
%       optional: display-flag
%output: cd= centroid distance function
%        centroid=[x y]
%        sample-points [x1 y1; x2 y2;..;xn yn]

%create edges from vertex-list P
last = length(P(:,1));
P=[P;[P(1,1) P(1,2)]];  %closed curve: connect endpoint to startpoint

i=1 : last;
vx = P(i+1,1)-P(i,1);
vy = P(i+1,2)-P(i,2);

%length
i=1:last;
l=sqrt(vx(i) .* vx(i) + vy(i) .* vy(i));

%cumulated length
cl=[];
for i=1:last;
   cl=[cl;sum(l(1:i))];
end


%distance between sample points along boundary
sampleDistance=cl(last)/n; 

%start with polygon's first vertex
cdx=P(1,1);
cdy=P(1,2);

%walk along boundary to find sample points
sampleIndex=2;
boundaryIndex=1;
while sampleIndex <= n
   while cl(boundaryIndex) < (sampleIndex-1) * sampleDistance
      boundaryIndex = boundaryIndex + 1;
   end
   %we might have moved too far, linear interpolation between the
   %last two boundary points results in the correct point
   toofar=cl(boundaryIndex) - (sampleIndex-1)*sampleDistance;
   backFactor=1- (toofar / l(boundaryIndex));
   
   cdx=[cdx; P(boundaryIndex,1) + (P(boundaryIndex+1,1)-P(boundaryIndex,1))* backFactor];   
   cdy=[cdy; P(boundaryIndex,2) + (P(boundaryIndex+1,2)-P(boundaryIndex,2))* backFactor];      
   
   %next sample
   sampleIndex=sampleIndex+1;
end

samples=[cdx cdy];
%now cdx,cdy contain sample points
%centroid
centroidX=sum(P(:,1))/last;
centroidY=sum(P(:,2))/last;
cdx=cdx - centroidX; %center
cdy=cdy - centroidY;

%distance
cd = sqrt(cdx .* cdx + cdy .* cdy);
centroid=[centroidX centroidY];

