# bwdistsc 快速距离场计算函数解析

## 声明与输入解析

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


## 计算距离场

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


### 计算一维距离场

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的循环结束了


### 计算三维距离场

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

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