bwdistsc 快速距離場計算函數解析

以前的《距離場計算：維度誘導法（dimension-induction）的基本原理》，已經簡述了快速距離場計算基本思路，bwdistsc函數是 Yuriy Mishchenko 對上述思路的實現。這裏將函數源碼作分析與解釋。介意結合 Yuriy Mishchenko 論文閱讀。

經測試，函數速度至關快，而且計算量與實際問題複雜度相關。若實際問題中物體較爲規則，在相同網格數量下，算法執行速度將比複雜拓撲結構物體相比有明顯提升。

文章爲簡明刪去了部分原有註釋，主要是權利聲明等。須要使用函數的話可到上文連接處Mathworks網站下載最新版

聲明與輸入解析

matlabfunction D=bwdistsc(bw,aspect)

% parse inputs
if(nargin<2 || isempty(aspect)) aspect=[1 1 1]; end

% determine geometry of the data
if(iscell(bw)) shape=[size(bw{1}),length(bw)]; else shape=size(bw); end

% correct this for 2D data
if(length(shape)==2) shape=[shape,1]; end
if(length(aspect)==2) aspect=[aspect,1]; end

% allocate internal memory
D=cell(1,shape(3)); for k=1:shape(3) D{k}=zeros(shape(1:2)); end

以上都是預處理。二維仍是三維，網格長寬比是否爲1，等等。

計算距離場

matlab%%%%%%%%%%%%% scan along XY %%%%%%%%%%%%%%%%
for k=1:shape(3)    %循環，按第三個座標（z）循環
    if(iscell(bw)) bwXY=bw{k}; else bwXY=bw(:,:,k); end
    %切片儲存在bw中

    % initialize arrays
    DXY=zeros(shape(1:2));
    D1=zeros(shape(1:2));

    % if can, use 2D bwdist from image processing toolbox    
    if(exist('bwdist') && aspect(1)==aspect(2))
        D1=aspect(1)^2*bwdist(bwXY).^2;
    else    % if not, use full XY-scan

計算一維距離場

切成一條一條的計算。一行一行，從第一行向最後一行掃略（從上到下），再從最後一行向第一行掃略（從下到上），從而生成了兩個記錄矩陣。這兩個矩陣中每個網格點記錄的值是當前列中，按當前方向掃略的狀況下，離本身最近的物體網格點。

下文中"on"-pixel指的是輸入bw中數值爲1的網格，即認爲是有物體的網格
雖然有不少「條」，可是都用了向量化計算，因此這裏沒有循環語句

matlab
% z的循環還在進行
%%%%%%%%%%%%%%% X-SCAN %%%%%%%%%%%%%%%        
% reference for nearest "on"-pixel in bw in x direction down

%  scan bottow-up (for all y), copy x-reference from previous row 
%  unless there is "on"-pixel in that point in current row, then 
%  that is the nearest pixel now
xlower=repmat(Inf,shape(1:2)); 

xlower(1,find(bwXY(1,:)))=1;    % fill in first row
for i=2:shape(1)
    xlower(i,:)=xlower(i-1,:);  % copy previous row
    xlower(i,find(bwXY(i,:)))=i;% unless there is pixel
end

% reference for nearest "on"-pixel in bw in x direction up
xupper=repmat(Inf,shape(1:2));

xupper(end,find(bwXY(end,:)))=shape(1);
for i=shape(1)-1:-1:1
    xupper(i,:)=xupper(i+1,:);
    xupper(i,find(bwXY(i,:)))=i;
end

% build (X,Y) for points for which distance needs to be calculated
idx=find(~bwXY); [x,y]=ind2sub(shape(1:2),idx);

% update distances as shortest to "on" pixels up/down in the above
% 將兩個矩陣合併起來
DXY(idx)=aspect(1)^2*min((x-xlower(idx)).^2,(x-xupper(idx)).^2);

計算二維距離場

matlab
        % z的循環還在繼續
        %%%%%%%%%%%%%%% Y-SCAN %%%%%%%%%%%%%%%
        % this will be the envelop of parabolas at different y
        D1=repmat(Inf,shape(1:2));

        p=shape(2);
        for i=1:shape(2)
            % some auxiliary datasets
            % 取出平面上的一列。從左向右按列掃描
            % d0中存放的是距離的平方
            d0=DXY(:,i);

            % selecting starting point for x:
            % * if parabolas are incremented in increasing order of y, 
            %   then all below-envelop intervals are necessarily right-
            %   open, which means starting point can always be chosen 
            %   at the right end of y-axis
            % * if starting point exists it should be below existing
            %   current envelop at the right end of y-axis
            dtmp=d0+aspect(2)^2*(p-i)^2;
                %均勻網格下，aspect的三個分類份量均爲1，能夠忽略。寫在這裏真影響理解……

            L=D1(:,p)>dtmp;    % 比較D1中的一列與dtmp的大小
                %一維數組L中儲存的是大小比較的結果量（0或1，認爲是bool）
            idx=find(L);            
            D1(idx,p)=dtmp(L);    % D1的最後一列被設置爲dtmp的一部分（取最小值）
            % 上面這幾句還能夠精簡


            % these will keep track along which X should 
            % keep updating distances            
            map_lower=L;
            idx_lower=idx;    %儲存了dtmp比D1小的那些位置

            % scan from starting points down in increments of 1
            % 從尾巴開始向前掃描，重複相似上面的比較
            % 可是隻關心lower位置，下面有解釋
            for ii=p-1:-1:1
                % new values for D
                dtmp=d0(idx_lower)+aspect(2)^2*(ii-i)^2;

                % these pixels are to be updated
                L=D1(idx_lower,ii)>dtmp;
                D1(idx_lower(L),ii)=dtmp(L);

                % other pixels are removed from scan
                map_lower(idx_lower)=L;                
                idx_lower=idx_lower(L);

                if(isempty(idx_lower)) break; end
            end
        end
    end
    D{k}=D1; 
end     %z的循環結束了

從尾巴開始向前掃描，重複相似上面的比較，可是只關心已經包含在lower中的位置。循環中不斷縮小lower區域，lower若爲空則能夠提早結束循環

緣由解釋

參考論文中引理1,拋物線比下包絡線低的位置只可能存在於一個區間，而不會斷斷續續的存在於多個區間。又因爲全部的拋物線二次項係數相等，因此這一區間一定爲無窮區間$(-\infty,x)$或者$(x,\infty)$或者$(-\infty,\infty)$（排除徹底重合的狀況）。若是是$(x,\infty)$的狀況，則下面這個從右向左掃描，找到包絡線與拋物線交點即再也不繼續的算法容易理解。

因爲上面的i循環是從左往右掃描的，新加入的拋物線要麼在包絡線右半部分上升段與之相交（上述區間爲$(x,\infty)$），或者整條都低於包絡線（上述區間爲$(-\infty,\infty)$）。$(-\infty,x)$的狀況不存在，由於它已經被掃描方法排除了。能夠畫圖觀察。

計算三維距離場

matlab%%%%%%%%%%%%% scan along Z %%%%%%%%%%%%%%%%
% 新建一個跟D同樣大的元胞數組D1，用Inf填充
D1=cell(size(D));
for k=1:shape(3) 
  D1{k}=repmat(Inf,shape(1:2)); 
end

% start building the envelope 
p=shape(3);
for k=1:shape(3)
    % if there are no objects in this slice, nothing to do
    if(isinf(D{k}(1,1)))    %有一個是Inf，即說明這一層裏面一個物體點都沒找到
      continue;
    end

下面的算法，跟由一維構建二維方式相同。因爲每次在一個維度上作文章，因此不管擴展到多少維度，須要計算的老是隻有那麼一小撮拋物線，因此算法基本同樣。須要的話都能寫出遞歸來。

matlab% selecting starting point for (x,y):
    % * if parabolas are incremented in increasing order of k, then all 
    %   intersections are necessarily at the right end of the envelop, 
    %   and so the starting point can be always chosen as the right end
    %   of the axis

    % check which points are valid starting points, & update the envelop
    dtmp=D{k}+aspect(3)^2*(p-k)^2;
    L=D1{p}>dtmp; 
    D1{p}(L)=dtmp(L);    

    % map_lower keeps track of which pixels can be yet updated with the 
    % new distance, i.e. all such XY that had been under the envelop for
    % all Deltak up to now, for Deltak<0
    map_lower=L;

    % these are maintained to keep fast track of whether map is empty
    idx_lower=find(map_lower);

    % scan away from the starting points in increments of -1
    for kk=p-1:-1:1
        % new values for D
        dtmp=D{k}(idx_lower)+aspect(3)^2*(kk-k)^2;

        % these pixels are to be updated
        L=D1{kk}(idx_lower)>dtmp;
        map_lower(idx_lower)=L;
        D1{kk}(idx_lower(L))=dtmp(L);

        % other pixels are removed from scan
        idx_lower=idx_lower(L);

        if(isempty(idx_lower)) break; end
    end
end

% prepare the answer
if(iscell(bw))
    D=cell(size(bw));
    for k=1:shape(3) D{k}=sqrt(D1{k}); end
else
    D=zeros(shape);
    for k=1:shape(3) D(:,:,k)=sqrt(D1{k}); end
end
end

計算完成