function generator = mpb ( ndim, n, generator, npp )
%*****************************************************************************80
%
%% MPB projects generators onto the boundary of the region.
%
% Discussion:
%
% The number NPP sets the number of subintervals into which we subdivide
% the boundary. It does NOT specify how many points will be pulled onto
% the boundary. The reason for this is that, after the first boundary
% subinterval has had a generator pulled into it, on every subsequent
% subinterval the nearest generator is likely to be the one in the
% previous subinterval! Unless an interior generator is closer than
% some small distance, this process will simply drag some unfortunate
% generator onto the boundary, and then around from interval to interval
% for a considerable time.
%
% The algorithm could be changed, if desired, so that points snapped
% to the boundary are guaranteed not to move, at least not twice in
% one application of this routine!
%
% Licensing:
%
% This code is distributed under the GNU LGPL license.
%
% Modified:
%
% 02 December 2004
%
% Author:
%
% Lili Ju
%
% Input:
%
% integer NDIM, the spatial dimension.
%
% integer N, the number of generators.
%
% real GENERATOR(NDIM,N), the coordinates of
% the generators. On output, some generators will have been moved.
%
% integer NPP, the number of subintervals into which the
% perimeter is divided.
%
% Output:
%
% real GENERATOR(NDIM,N), the updated coordinates of the generators.
%
dx = 1.0;
dy = 1.0;
%
% HH is the length of an individual segment of the perimeter of the region.
%
% U is set in such a way that on step I, it measures the distance from
% the lower left corner of the box to the midpoint of the I-th subinterval
% on the perimeter of the box.
%
hh = 2.0 * ( dx + dy ) / npp;
u = -0.5 * hh;
for i = 1 : npp
u = u + hh;
%
% The portion of the bottom perimeter from (0,0) to (1,0).
%
if ( u < dx )
sample(1:2,1) = [ u, 0.0 ]';
nearest = find_closest ( ndim, n, 1, sample, generator );
generator(2,nearest(1)) = 0.0;
%
% The portion of the right perimeter from (1,0) to (1,1).
%
elseif ( dx < u && u < dx + dy )
sample(1:2,1) = [ 1.0, ( u - dx ) ]';
nearest = find_closest ( ndim, n, 1, sample, generator );
generator(1,nearest(1)) = 1.0;
%
% The portion of the top perimeter from (1,1) to (0,1).
%
elseif ( dx + dy < u && u < dx + dy + dx )
sample(1:2,1) = [ 1.0 - ( u - dx - dy ), 1.0 ]';
nearest = find_closest ( ndim, n, 1, sample, generator );
generator(2,nearest(1)) = 1.0;
%
% The portion of the left perimeter from (0,1) to (0,0).
%
elseif ( dx + dy + dx < u )
sample(1:2,1) = [ 0.0, 1.0 - ( u - dx - dy - dx ) ]';
nearest = find_closest ( ndim, n, 1, sample, generator );
generator(1,nearest(1)) = 0.0;
end
end
return
end