function varargout = distributePointsOnPolygon( varargin )
% DISTRIBUTEPOINTSONPOLYGON this function uses faces and 2D/3D vertices as
% input. The vertices are redistributed over the polygon to decrease the
% variability in distance between the points.
% This is done by shifting the points iteratively over the original contour
% described by the polygon.
%
% Notes:
% - the polygon should be connected; i.e. the first and last point must be
% connected by faces.
% - the vertices may not be connected > 2 times.
% - the resulting polygon might be slightly smoohted, depending on the
% number of original points and the distance between those points; you can
% avoid this by using more input points or adding points at the longest faces
%
% Usage:
% [ v ] = distributePointsOnPolygon ( v )
% [ f, v ] = distributePointsOnPolygon ( f, v )
% [ f, v ] = distributePointsOnPolygon ( f, v, debug )
%
% use debug = 1 to see the progress.
% tweak the PARAMETRS on line 70-76 to speed up the process/ improve accuracy
%
%
% Demo:
% distributePointsOnPolygon()
%
%
% --- Author information
% Wouter Potters
% Academic Medical Center, Amsterdam, The Netherlands
% w.v.potters@amc.nl
% Date: 4-June-2013
% check input and output arguments
narginchk(0,3);
nargoutchk(0,2);
switch nargin
case 0
% demo time!!!
demoFunction();
return;
case 1 % assume faces [1:end] and [2:end 1]
v = varargin{1};
f = [ (1:size(v,1))' [2:size(v,1) 1]' ];
debug = 0;
case 2
f = varargin{1};
v = varargin{2};
debug = 0;
case 3
f = varargin{1};
v = varargin{2};
debug = varargin{3};
end
if debug == 1
disp('Debugging on; please note that the function takes longer when debugging is on')
end
if ( ( size(v,2) ~= 2 ) && ( size(v,2) ~= 3 ) ) || ( size(f,2) ~= 2 )
error('Vertices and faces should describe a 2D polygon (but may vertices may be 3D).')
end
v2 = v; f2 = f; % save original arguments
% add a point to connect the endpoints if required
if (f2(end,end) ~= f2(1,1)) && (f2(end,end) ~= f2(1,end)) && (f2(end,1) ~= f2(1,1)) && (f2(end,1) ~= f2(1,end))
f2(end+1,:) = [f2(end,end) f2(1,1)]; % add 1 connecting face
end
%%% PARAMETERS START
threshold = 0.10; % initial threshold
threshold_steps = 0.04; % steps for threshold
final_threshold = 0.90; % final threshold
move_scaling = 0.8; % this scales the movement of a point to its next position
maximum_iterations = 100;
%%% PARAMETERS END
continue_updates = 1; % while loop condition
iteration_number = 1;
while continue_updates % is aborted when dividing 2 distances
% calculate length of faces = distance between points
dist_f2 = (sqrt(sum((v2(f2(:,1),:) - v2(f2(:,2),:)).^2,2)));
% calculate distances to neighboring points
diff_v = zeros(length(v2),2);
row_in_f = zeros(length(v2),2);
for iv = 1:length(v2)
[row,~] = find(f2 == iv);
if length(row) ~= 2
error('error length r ~= 2')
end
row_in_f(iv,:) = row; clear temp;
diff_v(iv,:) = [dist_f2(row_in_f(iv,1)),dist_f2(row_in_f(iv,2))];
end
% Calculate relative difference between connected faces - force value
% between 0 and 1
divided_diff_v = diff_v(:,1) ./ diff_v(:,2);
divided_diff_v( divided_diff_v > 1 ) = 1 ./ divided_diff_v( divided_diff_v > 1 );
% define 'badly' distributed points based on threshold
change_pts = find( divided_diff_v < threshold );
% plot progress - show which points will change
if debug == 1
debugPlotFunction( v, v2, divided_diff_v, change_pts, threshold );
end
% calculate new position for current (ic) vertex
for ic = change_pts'
current_total_distance = sum( diff_v(ic,:) );
smallest_distance = min( diff_v(ic,:) );
biggest_distance = max( diff_v(ic,:) );
selection = row_in_f( ic, ...
diff_v(ic,:) == biggest_distance ); % select target positon based on largest face
id_of_target = f2( selection, ...
f2( selection , : ) ~= ic ); % get index of vertex from facelist
distance_to_go = ( (current_total_distance / 2) - smallest_distance ) / biggest_distance;
direction_to_go = ( v2( id_of_target, : ) - v2( ic, : ) );
vector = move_scaling * distance_to_go * direction_to_go ; % direction and size of movement
v2( ic, : ) = v2( ic, : ) + vector; % move v2(ic,:) according to vector
end
if threshold > final_threshold
if ~any(divided_diff_v < threshold)
continue_updates = 0; % abort loop
end
else
threshold = threshold + threshold_steps; % increase threshold
end
if iteration_number > maximum_iterations
error('Maximum iterations reached')
else
iteration_number = iteration_number + 1;
if debug
disp(['Now performing iteration #' num2str(iteration_number) ' - current threshold: ' num2str(threshold)])
end
end
end
if debug == 1
debugPlotFunction( v, v2, divided_diff_v, [], threshold );
end
% set outputs
switch nargout
case 1
varargout{1} = v2;
case 2
varargout{1} = f2;
varargout{2} = v2;
end
end % function distributePointsOnPolygon
%% ADDITIONAL FUNCTIONS
%- DEBUGPLOTFUNCTION
function fi = debugPlotFunction( v0, v, divided_diff_v, change_pts, threshold)
fi = figure(1); set(fi,'color','w')
sp1 = subplot(1,2,1); cla;
switch size(v,2)
case 2
plot(v0(:,1),v0(:,2),'g.'); hold on
plot(v(:,1),v(:,2),'ro')
case 3
plot3(v0(:,1),v0(:,2),v0(:,3),'g.')
plot3( v(:,1), v(:,2), v(:,3),'ro')
end
axis equal tight off
view(2)
hold on;
if ~isempty(change_pts)
switch size(v,2)
case 2
plot(v(change_pts,1),v(change_pts,2),'b+')
case 3
plot3(v(change_pts,1),v(change_pts,2),v(change_pts,3),'b+')
end
lg = legend(sp1,'original points','current points','point subject to change in this iteration');
set(lg,'Location','NorthOutside')
end
subplot(1,2,2); cla;
hist(divided_diff_v,10);
c = hist(divided_diff_v,10);
axis([0 1 0 1.5*max(c)])
hold on
plot([threshold threshold],[0 max(c)],'r--')
title('Histogram: length of 1st and 2st face divided')
drawnow; % update graph
end % function debugPlotFunction
%- DEMOFUNCTION
function demoFunction()
answer = questdlg('Do you wish to draw your own polygon points?','DEMOFUNCTION FOR DISTRIBUTEPOINTSONPOLYGON','Yes','No','No');
switch answer
case 'Yes'
fh = figure(1313); axis([0 10 0 10]);
title('Please make a closed polygon of some points')
[v] = impoly();
try v = v.getPosition; catch, error('Poly failed; try again.'); end
if length(v) < 5
error('Please select at least 5 points')
end
% now increase the number of points to avoid smoothing of the
% contour
v = [ v; v(1,:)]; % add first point before interpolation
v2 = zeros(length(v)*5,2);
v2(:,1) = interp1(1:length(v),v(:,1),linspace(1,length(v),length(v)*5));
v2(:,2) = interp1(1:length(v),v(:,2),linspace(1,length(v),length(v)*5));
v = v2(1:end-1,:); % remove last point after interpolation to avoid duplicate points
try close(fh); end % try to close the figure window.
% create faces belonging to the vertices
f = [(1:size(v,1))' [2:size(v,1) 1]' ];
case 'No'
v = [-21.2341 18.9143 -0.4923; -21.1714 19.3139 -0.4760
-21.1404 18.6266 -0.5040; -21.1404 19.3731 -0.4736
-21.1046 18.5146 -0.5086
-21.0156 18.1149 -0.5249
-21.0111 19.7136 -0.4597
-20.9628 17.7153 -0.5412
-20.8586 17.3156 -0.5575
-20.8479 20.1133 -0.4434
-20.7404 17.1215 -0.5654
-20.7404 20.2266 -0.4388
-20.5938 16.9159 -0.5738
-20.5716 20.5129 -0.4271
-20.3404 16.5325 -0.5894
-20.3404 20.7760 -0.4164
-20.3235 16.5163 -0.5901
-20.2054 20.9126 -0.4108
-19.9404 16.2182 -0.6022
-19.9404 21.1424 -0.4014
-19.8198 16.1166 -0.6063
-19.7167 21.3123 -0.3945
-19.5404 15.9387 -0.6136
-19.5404 21.4561 -0.3886
-19.2952 15.7169 -0.6226
-19.1404 15.5690 -0.6287
-19.1404 21.5608 -0.3844
-18.7894 15.3173 -0.6389
-18.7404 15.2796 -0.6405
-18.7404 21.5791 -0.3836
-18.3404 15.1024 -0.6477
-18.3404 21.6371 -0.3813
-18.0091 21.7119 -0.3782
-17.9404 15.0633 -0.6493
-17.9404 21.7330 -0.3774
-17.9078 21.7119 -0.3782
-17.5404 15.0953 -0.6480
-17.5404 21.5781 -0.3837
-17.1404 15.1217 -0.6469
-17.1404 21.5629 -0.3843
-16.7404 15.1836 -0.6444
-16.7404 21.5248 -0.3858
-16.5346 15.3173 -0.6389
-16.3404 15.4492 -0.6336
-16.3404 21.4559 -0.3887
-16.0155 21.3123 -0.3945
-15.9404 15.6387 -0.6258
-15.9404 21.2698 -0.3962
-15.8272 15.7169 -0.6226
-15.5404 15.9218 -0.6143
-15.5404 21.0013 -0.4072
-15.4566 20.9126 -0.4108
-15.3008 16.1166 -0.6063
-15.2483 20.5129 -0.4271
-15.1404 16.2758 -0.5999
-15.1404 20.3836 -0.4324
-14.9339 16.5163 -0.5901
-14.9029 20.1133 -0.4434
-14.7521 16.9159 -0.5738
-14.7404 16.9339 -0.5730
-14.7404 19.9332 -0.4507
-14.5494 19.7136 -0.4597
-14.5325 17.3156 -0.5575
-14.4762 17.7153 -0.5412
-14.4442 18.1149 -0.5249
-14.3404 18.4358 -0.5118
-14.3404 19.4301 -0.4713
-14.3124 18.5146 -0.5086
-14.2728 19.3139 -0.4760
-14.2667 18.9143 -0.4923];
f = [14 16; 12 14
20 22
27 30
32 33
35 36
38 40
42 45
46 48
51 52
54 56
61 62
69 70
68 66
66 65
64 63
63 60
59 57
55 53
50 49
44 43
43 41
39 37
31 29
29 28
26 25
23 21
19 17
13 11
11 9
3 1
2 4
1 2
5 3
6 5
8 6
9 8
15 13
17 15
21 19
25 23
28 26
34 31
37 34
41 39
47 44
49 47
53 50
57 55
60 59
65 64
70 68
67 69
62 67
58 61
56 58
52 54
48 51
45 46
40 42
36 38
33 35
30 32
24 27
22 24
18 20
16 18
10 12
4 7
7 10];
end
distributePointsOnPolygon(f,v,1);
end % function demoFunction
