数字图像处理MATLAB学习笔记（十）

数字图像处理MATLAB学习笔记（十）

Representation and Description 表示与描述

令 $S$ 表示一幅图像中像素的子集。如果 $S$ 中组成它们所有像素间都存在通路，就说两个像素 $p$ 和 $q$ 在 $S$ 中是连接着的。

对于 $S$ 中的任何像素 $p$ ，在 $S$ 中连接 $p$ 的像素的集合被称为连通分量。

如果仅有一个连通分量，$S$ 将被称为连通集。

在一幅图像中，如果 $R$ 是连通集，那么像素的子集 $R$ 被称为图像的区域。

边界是点的连接集合。如果边界上的点形成顺时针或逆时针序列，就称边界上的点为有序的。如果边界上的每个点恰好有两个值为1的相邻像素点，而且它们不是4邻接的，边界将被称为最低限度连接的。

内点被定义为区域内除边界外的任意位置的点。

1. Background

1.1 Function for Extracting Regions and Their Boundaries 提取区域及其边界的函数

计算二值图像中所有的连通分量(区域)：function bwlabel

[L, num] = bwlabel(f, conn)

conn指定了期望的连通性（默认为8连接），num是找到的连通分量数，L是标记矩阵，对每个连通分量分配唯一的1倒num的整数。连通性的值可以影响监测区域的数量。

function bwperim：这个函数返回二值图像g，其中仅包含f中所有区域的周界(边界)像素。

g = bwperim(f, conn)

conn默认为4连接。

function bwboundaries：提取二值图像f中所有区域的真实边界坐标

B = bwboundaries(f, conn, options)

conn相对于边界本身，默认值为8连接，参数options可以是’holes’和’noholes’，默认为’holes’，提取区域和孔洞的边界，也可以提取嵌套在区域内的区域边界；'noholes’只能得到区域或其子区域的边界。

输出B是P x 1的单元数组，P是物体数和孔洞数。单元数组中的每个单元都包含一个np x 2的矩阵，其中的行是边界像素的行和列的坐标，np是相应区域的边界像素数。

[B, L] = bwboundaries(...)

这种情况下，L是标记矩阵，使用不同整数来标记f的每个元素，背景像素标记为0，区域和孔洞数由max(L(😃)给出。

[B, L, NR, A] = bwboundaries(...)

这种情况下，NR为找到的孔洞数，A为逻辑稀疏矩阵。

function bound2im：以np x 2坐标数组形式给定边界b，产生在b中坐标处为1、背景为0、尺寸为MxN的二值图像g。

g = bound2im(b, M, N)

典型的M=size(f, 1)，N=size(f, 2)。

function image = bound2im(b, M, N)
%BOUND2IM Converts a boundaxy to an image.
%   IMAGE = BOUND2IM (b) converts b, an np-by-2 array containing the
%   integer coordinates of a boundary, into a binary image with 1s
%   in the locations of the coordinates in b and 0s el sewhere. The
%   height and width of the image are equal to the Mmin + H and Nmin
%   + W, where Mmin = min(b(:,1)) - 1, N = min(b(:,2)) - 1, and H
%   and W are the height and width of the boundary. In other words,
%   the image created is the smallest image that will encompass the
%   boundary while maintaining the its original coordinate values.
% 
%   IMAGE = BOUND2IM(b, M, N) places the boundary in a region of
%   size M-by-N. M and N must satisfy the following conditions:
% 
%       M >= max (b(:,1)) - min(b(:,1)) + 1
%       N >= max (b(:,2)) - min(b(:,2)) + 1
% 
%   Typically, M = size(f,1) and N = size(f,2), where f is the
%   image from which the boundary was extracted. In this way, the
%   coordinates of IMAGE and f are registered with respect to each
%   other.

% Check input
if size(b, 2) ~= 2
    error('The boundary must be of size np-by-2');
end

% Make sure the coordinates are integers.
b = round(b);

% Defaults
if nargin == 1
    Mmin = min(b(:, 1)) - 1;
    Nmin = min(b(:, 2)) - 1;
    H = max(b(:, 1)) - min(b(:, 1)) + 1; % Height of boundary
    W = max(b(:, 2)) - min(b(:, 2)) + 1; % Width of boundary
    M = H + Mmin;
    N = W + Nmin;
end

% Create the image.
image = false(M, N);
linearIndex = sub2ind([M, N], b(:, 1), b(:, 2));
image(linearIndex) = 1;

2. Representation 表示

2.1 Chain Codes 链码

链码Chain Codes通过指定长度与方向的直线段的连接序列来表示边界。一般都是建立在4/8连接之上的线段。

每条线段的方向都如图所示进行编码，这种链码被称为Freeman Chain Codes。

差分：以逆时针为基础，从后向前依次进行减法运算。以4-链码为例，遵循$1-0=1$、$0-1=-1=3$、$2-0=2$、$0-2=-2=2$、$3-0=3$、$0-3=-3=1$。

例如：初识链码为10103322，一次差分为3133030。

在MATLAB中，使用function fchcode计算保存在数组b中的排过序的np x 2个点边界的集合的Freeman Chain Codes：

c = fchcode(b, conn, dir)

其中，conn指明连接方式，默认为8；dir指明输出链码的方向，‘same’(默认)与b中点的方向一致，'reverse’则相反。

输出的字段含义如下：

• c.fcc = Freeman chain code (1 x np)

• c.diff = c.fcc 的一阶差分码(1 x np)

• c.mm = 最小值的整数部分(1 x np)

• c.diffmm = c.mm的一阶差分码(1 x np)

• c.x0y0 = 链码的起点坐标(1x2)

源码如下：

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.xOyO   = 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. If it does not, an error is issued. See table
%   below.
% 
%   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 prog rams.
% 
%   FREEMAN CHAIN CODE REPRESENTATION The table on the left shows
%   the 8-connected Freeman chain codes cor responding to allowed
%   deltax, deltay pairs. An 8-chain is converted to a 4-chain if
%   (1) 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. See Fig. 12.2 for an explanation of the di rectional 4-
%   and 8-codes.
% 
%       ----------------------------- ---------------
%       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 cor responding 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.

% 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 bwboundaries.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];

% Check the curve for out-of-order points or breaks.
% 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 de1 tay 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, organizedin 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);
end

% If 4-connectivity is specified, check that all components 
% of fcc are 0, 2, 4, or 6.
if conn == 4
    if isempty(find(fcc == 1 || fcc == 3 || fcc == 5 ...
            || fcc == 7, 1))
        fcc = fcc ./ 2;
    else
        error('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);

% Obtain code of the integer of minimum magnitude.
c.mm = minmag(fcc);

% 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 the
% code.
cr = fliplr(fcc);

% Next, obtain the new code values by traversing the code in the
% 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)
% Finds the integer of min imum 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 min imum 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 min imum magnitude.
[M, N] = size(A);
J = (1:M)';
D(J, 1) = 0; % Reserve memory space for loop.
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)
% 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 di fference is found by counting the number of direction
% changes (in a counter-clockwise direction)  that separate two
% adj acent 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

Example：Freeman Chain Code及其变体的实例

分别对原图像(第一张)进行9x9平均模版的平滑处理和阈值处理。然后取二值图像的边界，利用10%图像宽度的网格进行分离，得到点图，并连接序列。

2.2 使用最小周长多边形的多边形近似

数字边界能够用多边形以任意精度近似。对于闭合曲线，当多边形的顶点数目与边界点数目相同，并且每个顶点与边界点-致时，近似是精确的。多边形近似的目标是用尽可能少的顶点去表示给定边界的形状。

其中最为有效的算法是最小周长多边形(Minimum-Perimeter Polygon, MPP)算法。

MPP算法：用一组连接单元闭合一条边界。单元的大小决定于多边形近似的精度。

优势：多边形近似比链码、边界分段更具有抗噪声干扰的能力。对封闭曲线而言，当多边形的线段数与边界上点数相等时，多边形可以完全准确的表达边界。

限制：在实际应用中，多边形近似的目的是用最少的线段来表示边界，并且能够表达原边界的本质形状。

一些假设：围绕边界的单元集合称为细胞综合体。凸顶点记做白-W，凹顶点的镜像记做黑-B。

观点如下：

• 由简单连接的细胞联合体限制的MPP不是self-intersecting的。
• 每个MPP的凸顶点是W顶点，但并不是边界上的每个W顶点都是MPP的顶点。
• 每个MPP的镜像凹顶点是B顶点，但并不是边界上的每个B顶点都是MPP的顶点。
• 所有的B顶点都在MPP上或外边，并且所有的W顶点都在MPP上或里边。
• 最主要的是在顶点序列中，包含在细胞联合体中的最左边的顶点总是MPP的W顶点。

设定有三个点的方向： a = ( x a , y a ) a=(x_a,y_a) a=(x