在测试代码时发现，自己的matlab无法执行Freeman链码函数：
boundaries和fchcode函数都无法正常运行：
需要在自己的工作目录中添加如下函数：
boundaries fchcode minmag codediff
代码如下：
function B=boundaries(BW,conn,dir)
%BOUNDARIES Trace object boundaries.
%B=BOUNDARIES(BW) traces the exterior boundaries of objects in the binary
%image BW.B is a p_by_1 cell array,where p is the number of objects in the
%image.Each cell contains a Q_by_2 matrix,each row of which contains the
%row and column coordinates of a boundary pixel.Q is the number of boundary
%pixels for the corresponding object.Object boundaries are traced in the
%clockwise direction.
%
%B=BOUNDARIES(BW,CONN) specifies the connectivity to use when tracing
%boundaries.CONN may be either 8 or 4.The default value for CONN is 8.
%
%B=BOUNDARIES(BW,CONN,DIR) specifies the direction used for tracing
%boundaries.DIR should be either 'cw'(trace boundaries clockwise) or
%'ccw'(trace boundaries counterclockwise).If DIR is omitted BOUNDARIES
%traces in the clockwise direction.
if nargin<3 %nargin定义在用户自定义函数体内，nargin返回
%用来调用函数的变量的个数。
dir='cw';
end
if nargin<2
conn=8;
end
L=bwlabel(BW,conn); %返回值是矩阵L，大小同BW，包含了BW中的连通部分的标记
%The number of objects is the maximum value of L.Initialize the cell
%array(元包数组)B so that each cell initially contains a 0_by_2 matrix.
numObjects=max(L(:));%寻找L中所作标记中的最大值，这个最大值实际上就对应着L中
% 包含的最多的连通部分的数目。
if numObjects>0
B={zeros(0,2)};%元包数组中仅包含一个元素。
B=repmat(B,numObjects,1);%将B进行numObjects*1个复制构成新的B。
else
B={};
end
%Pad label matrix with zeros.This lets us write the boundary_following loop
%without worrying about going off the edge of the image.
Lp=padarray(L,[1 1],0,'both');
%Compute the linear indexing offsets to take us from a pixel to its
%neighbors.
M=size(Lp,1);%SIZE(X,1) returns the number of rows.
if conn==8
%Order is N NE E SE S SW W NW.
offsets=[-1,M-1,M,M+1,1,-M+1,-M,-M-1];
else
%Order is N E S W.
offsets=[-1,M,1,-M];
end
%next_search_direction_lut is a lookup table.Given the direction from pixel
%k to pixel k+1,what is the direction to start with when examining the
%neighborhood of pixel k+1?
if conn==8
next_search_direction_lut=[8 8 2 2 4 4 6 6];
else
next_search_direction_lut=[4 1 2 3];
end
%next_direction_lut is a lookup table.Given that we just looked at neighbor
%in a given direction,which neighbor do we look at next?
if conn==8
next_direction_lut=[2 3 4 5 6 7 8 1];
else
next_direction_lut=[2 3 4 1];
end
%Values used for marking the starting and boundary pixels.
START=-1;
BOUNDARY=-2;
%Initialize scratch space in which to record the boundary pixels as well as
%follow the boundary.
scratch=zeros(100,1);
%Find candidate starting locations for boundaries.
[rr,cc]=find((Lp(2:end-1,:)>0)&(Lp(1:end-2,:)==0));
rr=rr+1;
for k=1:length(rr)
r=rr(k);
c=cc(k);
if (Lp(r,c)>0)&(Lp(r-1,c)==0)&isempty(B{Lp(r,c)})
%We've found the start of the next boundary.Compute its linear
%offset,record which boundary it is,mark it,and initialize the
%counter for the number of boundary pixels.
idx=(c-1)*size(Lp,1)+r;
which=Lp(idx);
scratch(1)=idx;
Lp(idx)=START;
numpixels=1;
currentpixel=idx;
initial_departure_direction=[];
done=0;
next_search_direction=2;
while ~done
%Find the next boundary pixel.
direction=next_search_direction;
found_next_pixel=0;
for k=1:length(offsets)
neighbor=currentpixel+offsets(direction);
if Lp(neighbor)~=0
%Found the next boundary pixel.
if (Lp(currentpixel)==START)&...
isempty(initial_departure_direction)
%We are making the initial departure from the starting
%pixel.
initial_departure_direction=direction;
elseif (Lp(currentpixel)==START)&...
(initial_departure_direction==direction)
% We are about to retrace our path.
%That means we're done.
done=1;
found_next_pixel=1;
break;
end
%Take the next step along the boundary.
next_search_direction=...
next_search_direction_lut(direction);
found_next_pixel=1;
numpixels=numpixels+1;
if numpixels>size(scratch,1)
%Double the scratch space.
scratch(2*size(scratch,1))=0;
end
scratch(numpixels)=neighbor;
if Lp(neighbor)~=START
Lp(neighbor)=BOUNDARY;
end
currentpixel=neighbor;
break;
end
direction=next_direction_lut(direction);
end
if ~found_next_pixel
%If there is no next neighbor,the object must just have a
%single pixel.
numpixels=2;
scratch(2)=scratch(1);
done=1;
end
end
%Convert linear indices to row_column coordinates and save in the
%output cell array.
[row,col]=ind2sub(size(Lp),scratch(1:numpixels));
B{which}=[row-1,col-1];
end
end
if strcmp(dir,'ccw')
for k=1:length(B)
B{k}=B{k}(end:-1:1,:);
end
end
function c = fchcode(b, conn, dir)
%FCHCODE Computes the Freeman chain code of a boundary.
% C = FCHCODE(B) computes the 8-connected Freeman chain code of a
% set of 2-D coordinate pairs contained in B, an np-by-2 array. C
% is a structure with the following fields:
%
% c.fcc = Freeman chain code (1-by-np)
% c.diff = First difference of code c.fcc (1-by-np)
% c.mm = Integer of minimum magnitude from c.fcc (1-by-np)
% c.diffmm = First difference of code c.mm (1-by-np)
% c.x0y0 = Coordinates where the code starts (1-by-2)
%
% C = FCHCODE(B, CONN) produces the same outputs as above, but
% with the code connectivity specified in CONN. CONN can be 8 for
% an 8-connected chain code, or CONN can be 4 for a 4-connected
% chain code. Specifying CONN=4 is valid only if the input
% sequence, B, contains transitions with values 0, 2, 4, and 6,
% exclusively.
%
% C = FHCODE(B, CONN, DIR) produces the same outputs as above, but,
% in addition, the desired code direction is specified. Values for
% DIR can be:
%
% 'same' Same as the order of the sequence of points in b.
% This is the default.
%
% 'reverse' Outputs the code in the direction opposite to the
% direction of the points in B. The starting point
% for each DIR is the same.
%
% The elements of B are assumed to correspond to a 1-pixel-thick,
% fully-connected, closed boundary. B cannot contain duplicate
% coordinate pairs, except in the first and last positions, which
% is a common feature of boundary tracing programs.
% FREEMAN CHAIN CODE REPRESENTATION
% The table on the left shows the 8-connected Freeman chain codes
% corresponding to allowed deltax, deltay pairs. An 8-chain is
% converted to a 4-chain if (1) if conn = 4; and (2) only
% transitions 0, 2, 4, and 6 occur in the 8-code. Note that
% dividing 0, 2, 4, and 6 by 2 produce the 4-code.
%
% ----------------------- ----------------
% deltax | deltay | 8-code corresp 4-code
% ----------------------- ----------------
% 0 1 0 0
% -1 1 1
% -1 0 2 1
% -1 -1 3
% 0 -1 4 2
% 1 -1 5
% 1 0 6 3
% 1 1 7
% ----------------------- ----------------
%
% The formula z = 4*(deltax + 2) + (deltay + 2) gives the following
% sequence corresponding to rows 1-8 in the preceding table: z =
% 11,7,6,5,9,13,14,15. These values can be used as indices into the
% table, improving the speed of computing the chain code. The
% preceding formula is not unique, but it is based on the smallest
% integers (4 and 2) that are powers of 2.
% Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
% Digital Image Processing Using MATLAB, Prentice-Hall, 2004
% $Revision: 1.6 $ $Date: 2003/11/21 14:34:49 $
% Preliminaries.
if nargin == 1
dir = 'same';
conn = 8;
elseif nargin == 2
dir = 'same';
elseif nargin == 3
% Nothing to do here.
else
error('Incorrect number of inputs.')
end
[np, nc] = size(b);
if np < nc
error('B must be of size np-by-2.');
end
% Some boundary tracing programs, such as boundaries.m, output a
% sequence in which the coordinates of the first and last points are
% the same. If this is the case, eliminate the last point.
if isequal(b(1, :), b(np, :))
np = np - 1;
b = b(1:np, :);
end
% Build the code table using the single indices from the formula
% for z given above:
C(11)=0; C(7)=1; C(6)=2; C(5)=3; C(9)=4;
C(13)=5; C(14)=6; C(15)=7;
% End of Preliminaries.
% Begin processing.
x0 = b(1, 1);
y0 = b(1, 2);
c.x0y0 = [x0, y0];
% Make sure the coordinates are organized sequentially:
% Get the deltax and deltay between successive points in b. The
% last row of a is the first row of b.
a = circshift(b, [-1, 0]);
% DEL = a - b is an nr-by-2 matrix in which the rows contain the
% deltax and deltay between successive points in b. The two
% components in the kth row of matrix DEL are deltax and deltay
% between point (xk, yk) and (xk+1, yk+1). The last row of DEL
% contains the deltax and deltay between (xnr, ynr) and (x1, y1),
% (i.e., between the last and first points in b).
DEL = a - b;
% If the abs value of either (or both) components of a pair
% (deltax, deltay) is greater than 1, then by definition the curve
% is broken (or the points are out of order), and the program
% terminates.
if any(abs(DEL(:, 1)) > 1) | any(abs(DEL(:, 2)) > 1);
error('The input curve is broken or points are out of order.')
end
% Create a single index vector using the formula described above.
z = 4*(DEL(:, 1) + 2) + (DEL(:, 2) + 2);
% Use the index to map into the table. The following are
% the Freeman 8-chain codes, organized in a 1-by-np array.
fcc = C(z);
% Check if direction of code sequence needs to be reversed.
if strcmp(dir, 'reverse')
fcc = coderev(fcc); % See below for function coderev.
end
% If 4-connectivity is specified, check that all components
% of fcc are 0, 2, 4, or 6.
if conn == 4
val = find(fcc == 1 | fcc == 3 | fcc == 5 | fcc ==7 );
if isempty(val)
fcc = fcc./2;
else
warning('The specified 4-connected code cannot be satisfied.')
end
end
% Freeman chain code for structure output.
c.fcc = fcc;
% Obtain the first difference of fcc.
c.diff = codediff(fcc,conn); % See below for function codediff.
% Obtain code of the integer of minimum magnitude.
c.mm = minmag(fcc); % See below for function minmag.
% Obtain the first difference of fcc
c.diffmm = codediff(c.mm, conn);
function z = minmag(c)
%MINMAG Finds the integer of minimum magnitude in a chain code.
% Z = MINMAG(C) finds the integer of minimum magnitude in a given
% 4- or 8-connected Freeman chain code, C. The code is assumed to
% be a 1-by-np array.
% The integer of minimum magnitude starts with min(c), but there
% may be more than one such value. Find them all,
I = find(c == min(c));
% and shift each one left so that it starts with min(c).
J = 0;
A = zeros(length(I), length(c));
for k = I;
J = J + 1;
A(J, :) = circshift(c,[0 -(k-1)]);
end
% Matrix A contains all the possible candidates for the integer of
% minimum magnitude. Starting with the 2nd column, succesively find
% the minima in each column of A. The number of candidates decreases
% as the seach moves to the right on A. This is reflected in the
% elements of J. When length(J)=1, one candidate remains. This is
% the integer of minimum magnitude.
[M, N] = size(A);
J = (1:M)';
for k = 2:N
D(1:M, 1) = Inf;
D(J, 1) = A(J, k);
amin = min(A(J, k));
J = find(D(:, 1) == amin);
if length(J)==1
z = A(J, :);
return
end
end
function d = codediff(fcc, conn)
%CODEDIFF Computes the first difference of a chain code.
% D = CODEDIFF(FCC) computes the first difference of code, FCC. The
% code FCC is treated as a circular sequence, so the last element
% of D is the difference between the last and first elements of
% FCC. The input code is a 1-by-np vector.
%
% The first difference is found by counting the number of direction
% changes (in a counter-clockwise direction) that separate two
% adjacent elements of the code.
sr = circshift(fcc, [0, -1]); % Shift input left by 1 location.
delta = sr - fcc;
d = delta;
I = find(delta < 0);
type = conn;
switch type
case 4 % Code is 4-connected
d(I) = d(I) + 4;
case 8 % Code is 8-connected
d(I) = d(I) + 8;
end
添加完以上函数后，测试代码可以顺利运行：
close all;clear all;clc;
I = [1 1 1 1;1 1 0 1;0 1 0 1;0 1 1 1];
g = boundaries(I,4); %追踪4连接目标边界
c = fchcode(g{:},4); %求4方向freeman链码
c
output：
c =
x0y0: [1 1] %c.x0y0显示代码开始处的坐标（1×2）
fcc: [0 0 0 3 3 3 2 2 1 1 2 1] %c.fcc表示Freeman链码(1×n)，边界点集大小为 (n×2)
diff: [0 0 3 0 0 3 0 3 0 1 3 3] % c.diff 代码 c.fcc的一阶差分（1×n）
mm: [0 0 0 3 3 3 2 2 1 1 2 1] %c.mm表示最小幅度的整数
diffmm: [0 0 3 0 0 3 0 3 0 1 3 3] %c.diffmm表示代码c.mm的一阶差分（1×n）