Matlab code of chaincodes

Clear;
Close all;

%
B = imread('freeman.bmp ');
[M, N] = size (B );
% 8-connectivity, counterclockwise
B = boundaries (B, 4, 'cww ');
D = cellfun ('length', B );
[Max_d, K] = max (d );
B = B {1 };
BIM = bound2im (B, M, N, min (B (:, 1), min (B (:, 2 )));
Imwrite (BIM, '4link chaincode .bmp ');
% Hold on
% C = fchcode (B );
% Text (C. x0y0 (2), C. x0y0 (1), '/leftarrow', 'color', 'R', 'linewidth', 16 );
%
Gridsep = 50;
[S, Su] = bsubsamp (B, gridsep );
CN = connectpoly (S (:, 1), S (:, 2 ));
Bim2 = bound2im (CN, M, N, min (CN (:, 1), min (CN (:, 2 )));
Imwrite (bim2, ['grid sample', 'gridsep = ', num2str (gridsep), '.bmp']);
C = fchcode (Su, 4 );
C. diffmm

 

Function [s, Su] = bsubsamp (B, gridsep)
% Bsubsamp subsample a boundary.
% [S, Su] = bsubsamp (B, gridsep) subsamples the boundary B
% Assigning each of its points to the grid node to which it is
% Closest. The grid is specified by gridsep, which is
% Separation in pixels between the grid lines. For example, if
% Gridsep = 2, there are two pixels in between grid lines. So,
% Instance, the grid points in the first row wocould be ),
% (), (),..., And similarly in the Y direction. The value
% Of gridsep must be an even integer. the boundary is specified
% A set of coordinates in the form of an np-by-2 array. It is
% Assumed that the boundary is one pixel thick.
%
% Output S is the subsampled boundary. Output Su is normalized so
% That the grid separation is unity. This is useful for obtaining
% The Freeman chain code of the subsampled boundary.

% Copyright 2002-2004 R. C. gonzarez, R. E. Woods, & S. L. eddins
% Digital Image Processing Using MATLAB, Prentice-Hall, 2004
% $ Revision: 1.8 $ Date: 2004/11/04 20:17:59 $
 
% Check input.
[NP, NC] = size (B );
If NP <NC
Error ('B must be of size np-by-2 .');
End
If gridsep/2 ~ = Round (gridsep/2)
Error ('gridsep must be an even integer .')
End

% Some Boundary Tracing programs, such as boundaries. M, end
% The beginning, resulting in a sequence in which the coordinates
% Of the first and last points are the same. If this is the case
% In B, eliminate the last point.
If isequal (B (1, :), B (NP ,:))
NP = NP-1;
B = B (1: NP ,:);
End

% Find the max X and Y spanned by the boundary.
Xmax = max (B (:, 1 ));
Ymax = max (B (:, 2 ));

% Determine the Integral Number of grid lines with gridsep points in
% Between them that encompass the intervals [1, xmax], [1, Ymax].
Glx = Ceil (xmax + gridsep)/(gridsep + 1 ));
Ugly = Ceil (Ymax + gridsep)/(gridsep + 1 ));

% Form vectors of X and Y grid locations.
I = 1: Glx;
% Vector of Grid line locations intersecting x-axis.
X (I) = gridsep * I + (I-gridsep );

J = 1: ugly;
% Vector of Grid line locations intersecting Y-axis.
Y (j) = gridsep * j + (J-gridsep );

% Compute both components of the cityblock distance between each
% Element of B and all the grid-line intersections. Assign each
% Point to the grid location for which each comp of the cityblock
% Distance was <= gridsep/2. Note the use of meshgrid
% Optimize the code. Keep in mind that meshgrid requires that Columns
% Be listed first (see chapter 2 ).
Dist = gridsep/2;
[YG, XG] = meshgrid (Y, X );
Q = 1;
For k = 1: NP
[I, j] = find (ABS (XG-B (k, 1) <= DIST & ABS (YG-B (K, 2) <=...
Dist );
% IF point B (K, :) is equidistant from two or more grid intersections,
% Assign the point arbitrarily to the first one:
I = I (1 );
J = J (1 );
Ord = K; % to keep track of order of Input coordinates.
D1 (Q, = CAT (2, y (J), ORD );
D2 (Q, = CAT (2, X (I), ORD );
Q = q + 1;
End
 
% D is the set of points assigned to the new grid with line
% Separation of gridsep. Note that it is formed as D = (D2, D1)
% Compensate for the coordinate transposition inherent in using
% Meshgrid (see chapter 2 ).
D = CAT (2, D2 (:, 1), D1); % the second column of D1 & D2 is Ord.

% Sort the points using the values in ord, which is the last Col in
% D.
D = fliplr (d); % So the last column becomes first.
D = sortrows (d );
D = fliplr (d); % flip back.

% Eliminate duplicate rows in the first two components
% D to create the output. The CW or CCW order must be preserved.
S = D );
[S, m, n] = unique (S, 'rows ');

% Function unique sorts the data -- restore to original order
% By using the contents of M.
S = [s, m];
S = fliplr (s );
S = sortrows (s );
S = fliplr (s );
S = S );

% Scale to unit grid so that can use directly to obtain Freeman
% Chain code. The shape does not change.
% Su = round (S./gridsep) + 1;
SU (1, = Ceil (S (1, :)/gridsep) + 1;
For I = 1: length (S)-1;
If S (I + 1, 1)> S (I, 1)
DX = 1;
Elseif S (I + 1, 1) = S (I, 1)
DX = 0;
Else
DX =-1;
End
If S (I + 1, 2)> S (I, 2)
DY = 1;
Elseif S (I + 1, 2) = S (I, 2)
DY = 0;
Else
DY =-1;
End
SU (I + 1, = [SU (I, 1) + dx, SU (I, 2) + dy];
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
% 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 magn.pdf 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 abve,
% With the code connectivity specified in conn. Conn can be 8
% 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 abve,,
% In addition, the desired code direction is specified. Values
% 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
% 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, partition T 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 =
%,. These values can be used as indices into
% Table, improving the speed of computing the chain code.
% 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. gonzarez, 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 ('encrect 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
% 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.
% 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
% 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 abve.
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 magn.pdf.
C. Mm = minmag (FCC); % See below for function minmag.

% Obtain the first difference of FCC
C. diffmm = codediff (C. Mm, Conn );

% ------------------------------------------------------------------- %
Function Cr = coderev (FCC)
% Traverses the sequence of 8-connected Freeman chain code FCC in
% The opposite direction, changing the values of each code
% Segment. The starting point is not changed. FCC is a 1-by-np
% Array.

% Flip the array left to right. This redefines the Starting Point
% As the last point and reverses the order of "Travel" through
% Code.
Cr = fliplr (FCC );

% Next, obtain the new code values by traversing the code in
% Opposite direction. (0 becomes 4, 1 becomes 5,..., 5 becomes 1,
% 6 becomes 2, and 7 becomes 3 ).
Ind1 = find (0 <= CR & Cr <= 3 );
Ind2 = find (4 <= CR & Cr <= 7 );
Cr (ind1) = CR (ind1) + 4;
Cr (ind2) = CR (ind2)-4;

% ------------------------------------------------------------------- %
Function Z = minmag (c)
% Minmag finds the integer of minimum magnqueue in a chain code.
% Z = minmag (c) finds the integer of minimum magnen in a given
% 4-or 8-connected Freeman chain code, C. The code is assumed
% Be A 1-by-np array.

% The integer of minimum magn1_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
% Minimum magnloud. 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
% Elements of J. When length (j) = 1, one candidate remains. This is
% The integer of minimum magn.pdf.
[M, N] = size ();
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.
% Code FCC is treated as a circular sequence, so the last element
% Of D is the difference between the last and first elements
% FCC. The input code is a 1-by-np vector.
%
% The first difference is found by counting the number of ction
% 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