Farthest point sampling on 2d image

The principle of farthest point sampling is to randomly select a spot and then choose the point farthest from the point (the point with the highest value in D) and then continue iterating until you pick the number you want.

The main code is as follows:

%MAIN.M
clear options;
n = +;
[M,w] = Load_potential_map (' mountain ', n);

Npoints_list = Round (Linspace (20,200,6));% sample point number list
landmark = [];
Options.verb = 0;
ms = N;

CLF;
For I=1:length (npoints_list)
    nbr_landmarks = Npoints_list (i);
    Landmark = Perform_farthest_point_sampling (W, Landmark, Nbr_landmarks-size (landmark,2), options);%nbr_ Landmarks-size (landmark,2) minus the number of points that already exist
    %landmark the sampled points (including the original point and the new point)
    % compute the associated triangulation
    [D,z,q] = perform_fast_marching (W, landmark);
    % display sampling and distance function
    D = perform_histogram_equalization (D,linspace (0,1,n^2));% pull the value in D to [0,1] Within the range
    subplot (2,3,i);
    Hold on;
    Imageplot (D ');
    Plot (Landmark (1,:), Landmark (2,:), ' R. ', ' markersize ', MS);
    Title ([Num2str (Nbr_landmarks) ' points ']);
    Hold off;
    ColorMap Jet (n);
End

%PERFORM_FARTHEST_POINT_SAMPLING.M function [Points,d] = perform_farthest_point_sampling (W, points, npoints, options) 
% points is the point that has been sampled, npoints indicates the number of sample points that need to be added perform_farthest_point_sampling-samples points using farthest seeding strategy%
% points = perform_farthest_point_sampling (W, points, npoints);
% points can be is [] or can be a (2,npts) matrix of already computed% sampling locations.
% of Copyright (c) 2005 Gabriel Peyre options.null = 0;
If nargin<3 npoints = 1;
End If nargin<2 points = [];

End n = size (w,1);
Aniso = 0;
D = nb_dims (W);
    If d==4 aniso = 1; D = 2;
    % tensor field ElseIf d==5 aniso = 1; D = 3;
% tensor field end s = Size (W);

s = S (1:d);
% domain constraints (for shape meshing) L1 = getoptions (options, ' Constraint_map ', zeros (s) + Inf);
verb = getoptions (options, ' verb ', 1);

Mask = not (L1==-inf); If IsEmpty (points)% initialize farthest points at random points = round (rand (d,1) * (n-1)) +1;% Random one point D-dimensional coordinates% rEplace by farthest point [points,l] = perform_farthest_point_sampling (W, points, 1);% then point to points furthest distance Q = ones (s
    Ize (W)); Points = points (:, end);% takes the last point, that is, the point that is generated farthest from the initial random point npoints = npoints-1;% The number of points to generate minus 1 else% initial distance map [ L,Q] = My_eval_distance (W, points, options);% if some sample points are already present, you can calculate the distance by perform_fast_marching, points is the initial point (the point with a distance of 0)% L = min (
Zeros (s) + INF, L1);
% Q = zeros (s);
    End for i=1:npoints if npoints>5 && verb==1 ProgressBar (i,npoints);
    End options.nb_iter_max = INF; Options. Tmax = INF;
    % sum (size (W));
    % [D,s] = perform_fast_marching (W, points, options);
    Options.constraint_map = L;
    pts = Points;
    If not (aniso) pts = pts (:, end); Why does the% take only the most points? Because the distance in front is good, stored in L end D = My_eval_distance (W, pts, options);
    Dold = D; D = min (d,l); % known distance map to Lanmarks L = min (D,L1); % Cropp with and constraints if not (IsEmpty (Q))% update Voronoi Q (dold==d) = size (points,2);
    End% Remove Away data D (D==inf) = 0; If IsEmpty (Q)% compute farthest points [tmp,i] = max (D (:));% find the farthest point [a,b,c] = ind2sub (Size (W), I (1
    ));
        Else% compute farthest Steiner point [pts,faces] = compute_saddle_points (q,d,mask); A = PTS (n); B = pts (2,1);
        c = [];% 1th column, the maximum value for distance d if d==3 C = pts (3,1);
    End End If d==2% 2D points = [points,[a;b]];
    Else% 3D points = [points,[a;b;c]]; End End%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [d,q] = My_eval_distance (W, x, options)% gives the point weight matrix W,
The initial point x (the point with a distance of 0) calculates the distance of the points% D is distance% Q is Voronoi segmentation options.null = 0;
n = size (w,1);

D = nb_dims (W);
        If STD (w (:)) <eps% is the same as the W weight value Euclidean distance if size (x,2) >1 D = Zeros (n) +inf; Q =Zeros (n); For I=1:size (x,2) dold = D;
            Qold = Q;
            D = min (Dold, My_eval_distance (w,x (:, i)));
            % Update Voronoi segmentation Q (:) = i;
        Q (D==dold) = Qold (D==dold);
    End return;
        End If d==2 [y,x] = Meshgrid (1:n,1:n);
    D = 1/w (1) * sqrt ((x-x (1)). ^2 + (Y-x (2)). ^2);
        else [x, y, z] = Ndgrid (1:n,1:n,1:n);
    D = 1/w (1) * sqrt ((x-x (1)). ^2 + (Y-x (2)). ^2 + (Z-x (3)). ^2);
End Q = d*0+1;
else [d,s,q] = perform_fast_marching (W, X, Options); End

%PERFORM_FAST_MARCHING.M function [D,s,q] = perform_fast_marching (W, start_points, options)% Perform_fast_marching-la
Unch The Fast marching algorithm, in 2D or 3D.  
% [D,s,q] = perform_fast_marching (w, start_points, options)% W is a (n,n) (for A, d=2) or (N,n,n) (for the d=3) % weight matrix.
The geodesics would follow regions where W is large.
% W must be > 0.
% ' Start_points ' is a D x k array, start_points (:, i) are the ith starting point.
% D is the distance function to the set of starting points. % S is the final state of the points:-1 for dead (ie the distance% have been computed), 0 for open (ie the Dista NCE is only a temporary% value) and 1 for far (ie point not already computed).  Distance function% for remote points is Inf. (Note for the far point, 1 is the state, INF is the distance)% (according to the book,-1 for the known point, 0 for the trial point, 1 for the far points)% Q is the index of the closest point.
Q is the set to 0 for far points. 
% Q provide a Voronoi decomposition of the domain. %% OptIonal:%-can provide special conditions for stop in options:% ' options.end_points ': stop when these point
S was reached% ' Options.nb_iter_max ': stop when a given number of iterations is% reached. %-you can provide a heuristic in options.heuristic (typically. Try to guess the distance percent that remains fro
M a given node to a given target).       % this is an array of same size as W.%-can provide a map l=options.constraint_map that reduce the set of Explored points. Only points with current distance smaller than L% would be expanded.
Set some entries of L To-inf to avoid any% exploration of these points.
%-Options.values Set the initial distance value for starting points percent (default value is 0).
%% See also:perform_fast_marching_3d.

% of Copyright (c) Gabriel Peyre options.null = 0;
end_points = getoptions (options, ' end_points ', []);
verbose = getoptions (options, ' verbose ', 1); Nb_iter_max = geToptions (Options, ' Nb_iter_max ', INF);
Values = getoptions (options, ' values ', []);
L = getoptions (options, ' constraint_map ', []);
H = getoptions (options, ' heuristic ', []);

DMax = getoptions (options, ' DMax ', INF);

D = nb_dims (W); if (d==4 && size (w,3) ==2 && size (w,4) ==2) | | (d==4 && size (w,4) ==6) | |  (d==5 && size (w,4) ==3 && size (w,5) ==3)% anisotropic fast Marching if d==4 && size (w,3) ==2
        && size (w,4) ==2% 2D vector field, 3D field W1 = zeros (Size (w,1), size (w,2), 3, 3); 
        W1 (:,:, 1:2,1:2) = W;
        W1 (:,:, 3, 3) = 1;
        W = Reshape (W1, [Size (w,1) size (w,2), 1 3 3]); % convert to correct size W = Cat (4, W (:,:,:, Max), W (:,:,:, Max), W (:,:,:, 1,3), W (:,:,:, 2,2), W (:,:,:, 2,3), w (;        
    , 3, 3)); End If d==5% convert to correct size W = Cat (4, W (:,:,:, Max), W (:,:,:,), W (:,:,:, 1,3), W (:,:,:, 2,2
    ), W (:,:,:, 2,3), W (:,:,:, 3, 3)); End If SiZe (start_points,1) ==2 start_points (end+1,:) = 1;
    End If Size (start_points,1) ~=3 error (' start_points should be of size (3,n) ');
    End% Padd To avoid boundary problem W = Cat (1, W (1,:,:,:), W, W (End,:,:,:));
    W = Cat (2, W (:, 1,:,:), W, W (:, End,:,:));
    
W = Cat (3, W (:,:, 1,:), W, W (:,:, end,:));
 % if IsEmpty (l) L = ones (size (w,1), size (w,2), size (w,3));
    % End If dmax==inf dmax = 1E15;
    End% start_points = start_points-1;
    Alpha = 0;

    [D,q] = Perform_front_propagation_anisotropic (W, L, Alpha, Start_points,dmax);
    % Remove boundary problems d = d (2:end-1,2:end-1,2:end-1);
    Q = q (2:end-1,2:end-1,2:end-1);
    S = [];
    D (D>1E20) = INF;
Return
End If d~=2 && d~=3 error (' Works only in 2D and 3D. ');
End If Size (start_points,1) ~=d error (' Start_points should is (d,k) dimensional with k=2 or 3. ');
End L (l==-inf) =-1e9;
L (L==inf) =1e9; Nb_iter_max = min (Nb_iter_max, 1.2*max (Size (W)) ^d); % Use fast c-coded version if possible if d==2 if exist (' perform_front_propagation_2d ') ~=0 [d,s,q] = Perform_f
        Ront_propagation_2d (W,start_points-1,end_points-1,nb_iter_max, H, L, values); % of the VONOROI classification principle, assuming that the initial sample point has K, then the K sample point as the center of the K Cell, and then the rest of the point distance from which the sample point close to which cell% with the kind of plane cut out of the cell although the process is not the same, But the principle is the same. Each sample point will have a cell else [d,s] = Perform_front_propagation_2d_slow (w,start_points,end_points,nb_iter_
        Max, H);
    Q = [];

End ElseIf d==3 [d,s,q] = perform_front_propagation_3d (W,start_points-1,end_points-1,nb_iter_max, H, L, values);

End Q = q+1;
% replace C ' Inf ' Value (1e9) by Matlab Inf value. D (D>1e8) = INF;

The results of the operation are as follows:

Blue indicates a distance of 0, and red indicates a distance of 1.

Finally, the method is in common with medial axis:

1. The farthest point will definitely fall on the middle axis.

Proof: The furthest point is the distance of at least two points is equal, then the distance is farthest, then it must satisfy the condition that the distance is equal, that it will fall on the middle axis

2. It has a relationship with Power diagram: Powder diagram is inserted into the ball, equal to the area where the weight corresponding ball is set to 0. The smaller the weight value, the stronger the repulsion, the greater the attraction. If you set the area of the ball to 0.5, the resulting middle axis may be an arc, not a straight line. And the arc is relatively close to the value of the ball.

MATLAB Full source code