function voronoi_mountains ( p )
%*****************************************************************************80
%
%% voronoi_mountains() displays a surface plot of the Voronoi distance function.
%
% Discussion:
%
% For each point in the picture, we have computed the distance to the
% nearest Voronoi generator. This value is zero at a Voronoi generator,
% of course, and nearby, it rises in a circular cone. When two generators
% are neighbors, the boundary between them is a straight line segment;
% in a 3d plot, you see a ridge that is "straight" but whose height
% is parabolic.
%
% If you run this example, be sure to try out MATLAB's 3D graphing features,
% including the 3D rotation option and the zoom.
%
% Thanks to Jonathan Hadida for pointing out how to vectorize the
% distance calculation.
%
% Licensing:
%
% This code is distributed under the MIT license.
%
% Modified:
%
% 03 March 2011
%
% Author:
%
% John Burkardt
%
% Reference:
%
% Per-Olof Persson and Gilbert Strang,
% A Simple Mesh Generator in MATLAB,
% SIAM Review,
% Volume 46, Number 2, June 2004, pages 329-345.
%
% Input:
%
% real P(NP,2), the coordinates of a set of Voronoi generators.
%
% Local:
%
% integer NS, the number of grid points in the X and Y
% directions. Making this number bigger may give greater detail in the
% plot.
%
fprintf ( 1, '\n' );
fprintf ( 1, 'VORONOI_MOUNTAINS:\n' );
fprintf ( 1, ' MATLAB/Octave version %s\n', version ( ) ) ;
fprintf ( 1, ' Display a plot of a Voronoi region as a 3D surface.\n' );
fprintf ( 1, ' The Z coordinate in this plot is the distance of any point\n' );
fprintf ( 1, ' (X,Y) to the nearest Voronoi generator.\n' );
ns = 401;
xmin = min ( p(:,1) );
xmax = max ( p(:,1) );
xrange = xmax - xmin;
xmin = xmin - 0.30 * xrange;
xmax = xmax + 0.30 * xrange;
xdel = ( xmax - xmin ) / ( ns - 1 );
ymin = min ( p(:,2) );
ymax = max ( p(:,2) );
yrange = ymax - ymin;
ymin = ymin - 0.30 * yrange;
ymax = ymax + 0.30 * yrange;
ydel = ( ymax - ymin ) / ( ns - 1 );
%
% Set up NS sample points.
% Unfortunately, XS and YS are matrices, and we will shortly need them
% to be column vectors...and then matrices again!
%
[ xs, ys ] = meshgrid ( xmin : xdel : xmax, ymin : ydel : ymax );
%
% Compute the Delaunay triangulation of the data points,
% and then use the Delaunay search method.
%
% MATLAB's "delaunay" and "dsearch" functions are going to be retired soon.
%
if ( false )
t = delaunay ( p(:,1), p(:,2) );
k = dsearch ( p(:,1), p(:,2), t, xs, ys );
else
t = DelaunayTri ( p(:,1), p(:,2) );
xs = reshape ( xs, ns * ns, 1 );
ys = reshape ( ys, ns * ns, 1 );
k = nearestNeighbor ( t, xs, ys );
end
%
% Evaluate the distance.
% I am sure you can do this in a vectorized way, but I gave up trying.
%
if ( false )
for j = 1 : ns
for i = 1 : ns
zs(i,j) = sqrt ( ( xs(i,j) - p(k(i,j),1) ).^2 ...
+ ( ys(i,j) - p(k(i,j),2) ).^2 );
end
end
%
% Jonathan Hadida pointed out how to vectorize the above loop.
%
else
zs = sqrt ( ( xs - reshape ( p(k,1), ns * ns, 1 ) ).^2 ...
+ ( ys - reshape ( p(k,2), ns * ns, 1 ) ).^2 );
end
%
% Draw the 3D surface.
%
clf ( );
xs = reshape ( xs, ns, ns );
ys = reshape ( ys, ns, ns );
zs = reshape ( zs, ns, ns );
h = surf ( xs, ys, zs, 'FaceColor', 'Interp', 'EdgeColor', 'None' );
%
% Choose a color map and the number of levels.
% An interesting and different choice is "prism ( 60 )".
%
% colormap ( cool ( 60 ) );
colormap ( prism ( 60 ) );
hold ( 'on' );
%
% Add little dots to the plot indicating the location of the generators.
%
z = zeros ( size(p,1), 1 ) + 0.04;
scatter3 ( p(:,1), p(:,2), z, 100, 'w', 'filled' );
view ( 2 );
axis ( 'equal' );
axis ( 'off' );
hold ( 'off' );
return
end