Farthest Point Sampling on 2d image

Farthest Point Sampling的原理是,先隨機選一個點,然後呢選擇離這個點距離最遠的點(D中值最大的點)加入起點,然後繼續迭代,直到選出需要的個數為止

其主要代碼如下:

%main.mclear options;n = 400;[M,W] = load_potential_map('mountain', n);npoints_list = round(linspace(20,200,6));%採樣點個數列表landmark = [];options.verb = 0;ms = 15;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) 減去已經存在的點數    %landmark為已採樣的點(包括原來的點和新增的點)    % 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));%把D中的值拉到[0,1]範圍內    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(256);end

%perform_farthest_point_sampling.mfunction [points,D] = perform_farthest_point_sampling( W, points, npoints, options )% points為已經採樣了的點,npoints表示需要加入採樣點的個數% perform_farthest_point_sampling - samples points using farthest seeding strategy%% points = perform_farthest_point_sampling( W, points, npoints );%%   points can be [] or can be a (2,npts) matrix of already computed %       sampling locations.%   %   Copyright (c) 2005 Gabriel Peyreoptions.null = 0;if nargin<3    npoints = 1;endif nargin<2    points = [];endn = size(W,1);aniso = 0;d = nb_dims(W);if d==4    aniso = 1;    d = 2; % tensor fieldelseif d==5    aniso = 1;    d = 3; % tensor fieldends = 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;%隨機一個點的d維座標    % replace by farthest point    [points,L] = perform_farthest_point_sampling( W, points, 1 );%然後選點到points最遠的距離    Q = ones(size(W));    points = points(:,end);%取最後一個點,即就是產生的離初始隨機點最遠的那個點    npoints = npoints-1;%需要產生的點數減1else    % initial distance map    [L,Q] = my_eval_distance(W, points, options);%如果初始已存在一些採樣點,則可以通過perform_fast_marching算距離了, points為初始點(距離為0的點)%    L = min(zeros(s) + Inf, L1);%    Q = zeros(s);endfor 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);%為何只取最一個點?因為前面的距離都算好了,儲存在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 other 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(:));%找距離最遠的點        [a,b,c] = ind2sub(size(W),I(1));    else        % compute farthest steiner point        [pts,faces] = compute_saddle_points(Q,D,mask);        a = pts(1,1); b = pts(2,1); c = [];%第1列,為距離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)%給點權值矩陣W, 初始點x(距離為0的點),則算各點的距離% D is distance% Q is voronoi segmentationoptions.null = 0;n = size(W,1);d = nb_dims(W);if std(W(:))<eps%即W權值裡面值一樣    % 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.mfunction [D,S,Q] = perform_fast_marching(W, start_points, options)% perform_fast_marching - launch the Fast Marching algorithm, in 2D or 3D.%%   [D,S,Q] = perform_fast_marching(W, start_points, options)%%   W is an (n,n) (for 2D, d=2) or (n,n,n) (for 3D, d=3) %       weight matrix. The geodesics will follow regions where W is large.%       W must be > 0.%   'start_points' is a d x k array, start_points(:,i) is 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%       has been computed), 0 for open (ie the distance is only a temporary%       value), 1 for far (ie point not already computed). Distance function%       for far points is Inf.(注意對於far來說,1是狀態,Inf是距離)%       (按照書上的說法,-1為known的點,0為trial點,1為far點)%   Q is the index of the closest point. Q is set to 0 for far points.%       Q provide a Voronoi decomposition of the domain. %%   Optional:%   - You can provide special conditions for stop in options :%       'options.end_points' : stop when these points are reached%       'options.nb_iter_max' : stop when a given number of iterations is%          reached.%   - You can provide an heuristic in options.heuristic (typically that try to guess the distance%       that remains from a given node to a given target).%       This is an array of same size as W.%   - You can provide a map L=options.constraint_map that reduce the set of%       explored points. Only points with current distance smaller than L%       will 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%   (default value is 0).%%   See also: perform_fast_marching_3d.%%   Copyright (c) 2007 Gabriel Peyreoptions.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(:,:,:,1,1), W(:,:,:,1,2), W(:,:,:,1,3), W(:,:,:,2,2), W(:,:,:,2,3), W(:,:,:,3,3) );            end    if d==5        % convert to correct size        W = cat(4, W(:,:,:,1,1), W(:,:,:,1,2), 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;endif d~=2 && d~=3    error('Works only in 2D and 3D.');endif size(start_points,1)~=d    error('start_points should be (d,k) dimensional with k=2 or 3.');endL(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 possibleif d==2    if exist('perform_front_propagation_2d')~=0        [D,S,Q] = perform_front_propagation_2d(W,start_points-1,end_points-1,nb_iter_max, H, L, values);        %講下vonoroi的分類原理, 假設初始sample點有k個,那麼就把這k個sample點作為k個cell的中心,然後將剩下的點距離哪個sample點近就把歸於哪個cell裡        %跟那種用平面切出來的cell雖然過程不一樣,但是原理是一樣的.每一個sample點會擁有一個cell    else        [D,S] = perform_front_propagation_2d_slow(W,start_points,end_points,nb_iter_max, H);        Q = [];    endelseif d==3    [D,S,Q] = perform_front_propagation_3d(W,start_points-1,end_points-1,nb_iter_max, H, L, values);endQ = Q+1;% replace C 'Inf' value (1e9) by Matlab Inf value.D(D>1e8) = Inf;

運行結果如下:

藍色表示距離為0, 紅色表示距離為1.

最後講下該方法與medial axis的共同之處:

1. 最遠點一定會落在中軸上面

證明: 最遠點是指至少到兩個點的距離相等,則此距離最遠,那麼它肯定滿足距離相等這一條件,即它一定會落在中軸上面

2.它與power diagram的關係為:powder diagram插進球後, 相等於將把weight對應的球的地區設為0.  weight值越小,排斥力越強, 越大,吸引力越強.如果把整個球的地區設為0.5,那麼產生的中軸可能就是弧形,而不是直線.而且該弧開是比較靠近值大的球.

matlab完整原始碼