实验目的

通过编程，实现将一幅彩色图像分割为若干个同质区域，即采用K-Means聚类算法来将像素分组从而实现图像分割。在实验中，要分别基于颜色特征和纹理特征实现图像分割，并通过分析比较两种视觉特征在图像分割中的性能。

实验设计

本实验设计4个功能函数

生成隶属度矩阵

函数说明

给定一个h * w *d 的矩阵featIm，其中h 和w 为原始图像的高度和宽度，d 表示图像中每一个像素点所提取的特征向量的维数。给定一个有k 个聚类中心点的矩阵meanFeatures(矩阵维度是k*d)，其中每个中心点都是一个d 维的行向量（矩阵的一行），将输入图片中的每个像素映射到其所归属的k-means 某一个中心点去。函数的返回值定义为labelim，labelim 是一个h*w 的整数矩阵，用来表明每一个像素所属的聚类中心标号(1…k)。函数形式：

function [labelIm] = quantizeFeats (featIm, meanFeats)

实现思路

featIm 是一个 d 维矩阵，meanFeats 里面有 k 个中心点，每个中心点有 d 列，利用 kmeans 的思想，依次计算 featIm 中的每一个像素点到 k 个中心点的距离，选取最近的中心点的行号作为该点的标记值进行返回，返回的 labelIm 正好是 h*w 的矩阵。

代码展示

function labelIm = quantizeFeats (featIm, meanFeats)

height = size(featIm, 1);
width = size(featIm, 2);
dimension = size(featIm, 3);

% 计算原始图像 至 meanFeats 的距离，生成隶属度矩阵 labelIm
labelIm = ones(height, width);
for h=1:height
    for w=1:width
       minDis = -1;
       for k=1:size(meanFeats, 1)
           dis = 0;
           for d=1:dimension
               dis = dis + power(featIm(h, w, d) - meanFeats(k, d), 2);
           end
           dis = sqrt(dis);
           if minDis < 0 || dis < minDis
               labelIm(h, w) = k;
               minDis = sqrt(dis);
           end
       end
    end 
end

return

% 给定一个h * w *d 的矩阵featim，其中h 和w 为原始图像的高度和宽度，d 表示图像中每一个像素点所提取的特征向量的维数。
% 给定一个有k 个聚类中心点的矩阵meanFeatures(矩阵维度是k*d)，其中每个中心点都是一个d 维的行向量（矩阵的一行），
% 将输入图片中的每个像素映射到其所归属的k-means 某一个中心点去。
% 函数的返回值定义为labelim，labelim 是一个h*w 的整数矩阵，用来表明每一个像素所属的聚类中心标号(1...k)。

生成纹理基元编码集

函数说明

给定一个长度为n 且包含n 个灰度值图像的元胞数组imStack，以及滤波器组bank，基于所有n 个图片的过滤响应样本计算一个纹理基元编码集（如一组量化的滤波器组响应）。注意，元胞数组的特点是数组中的每个元素可以存储不同大小的矩阵，所以允许每张图片有不同的宽和高。其中bank 是一个包含d 个滤波器的m*m*d 矩阵，每个滤波器的大小为m*m，textons 是一个k*d 的矩阵，其中每一行代表一个纹理特征，如一个量化滤波器组的响应。函数形式：

function [textons] = createTextons(imStack, bank, k)

实现思路

imStack 是元胞数组，里面承载的是n 个灰度值图像，此处使用的滤波器组bank是49*49*38的矩阵，它包含38 个总过滤器，其中每个过滤器为49 x 49。首先将每一个灰度值图像使用滤波器组过滤，得到一个r1*c1*38的矩阵，其次把所有图像的过滤响应样本都合并为一个 row*38 的矩阵，再利用kmeans算法将其聚为 k 类，得到的聚类中心就可以作为一个纹理基元编码集。

代码展示

function [textons] = createTextons(imStack, bank, k)

[row, col] = size(imStack);
bankNum = size(bank, 3);

% 用滤波器过滤元胞数组中的n个灰度图，得到n个 bankNum 维的过滤响应样本
% 再将这n个过滤响应样本组合成一个大的样本集 textonsData
textonsData = [];
for i=1:row
    for j=1:col
        im = imStack{i, j};
        im = im2double(im);
        responses = zeros(size(im, 1), size(im, 2), bankNum);
        for r=1:bankNum
            responses(:,:,r)=conv2(im,double(bank(:,:,r)),'same');
        end
        X = reshape(responses, size(responses,1)* size(responses,2), bankNum);
        Xrow = size(X, 1);
        ranX = X(randperm(Xrow, ceil(Xrow/1000)),:);
        textonsData = [textonsData; ranX];
    end
end

% 使用 kmeans 方法聚合出含有k个纹理基元的纹理编码集
[~, textons] = kmeans(textonsData, k);
return;

% 给定一个长度为n 且包含n 个灰度值图像的元胞数组imStack，以及滤波器组bank，
% 基于所有n 个图片的过滤响应样本计算一个纹理基元编码集（如一组量化的滤波器组响应）。
% 其中bank 是一个包含d 个滤波器的m*m*d 矩阵，每个滤波器的大小为m*m，textons 是一个k*d 的矩阵，
% 其中每一行代表一个纹理特征，如一个量化滤波器组的响应。

构建纹理柱状图

函数说明

给定一张灰度图像，一个滤波器组，一份纹理编码集，构建纹理柱状图。对于每个像素，基于每个纹理在其邻近范围（定义在固定大小的winSize 内的局部窗口）内出现的频率。函数形式：

function [featIm] = extractTextonHists(origIm, bank, textons, winSize)

实现思路

将原始灰度图像origIm，使用滤波器组进行过滤，得到一个r*c*38 的矩阵，计算每个像素点到纹理基元编码集的隶属度矩阵feattexton。然后统计像素点的邻域窗口内，每个纹理基元出现的频率，可以用一个向量表示，设textons 中有 k 个纹理基元，那么这个向量就有k个数据，每个数据即为相应基元出现的频率，再将这个向量作为 featIm 中对应像素点的数据元素返回，就得到了直方图矩阵，即 featIm 为r*c*k 的矩阵。

代码展示

function [featIm] = extractTextonHists(origIm, bank, textons, winSize)

origIm = im2double(origIm);
[row, col] = size(origIm);

% 图像过滤，生成滤波器组响应 responses
bankNum = size(bank, 3);
responses = zeros(row, col, bankNum);
for r=1:bankNum
    responses(:,:,r)=conv2(origIm,double(bank(:,:,r)),'same');
end

% 计算 滤波器组响应 与 纹理基元编码集 的距离，并生成隶属度矩阵 feattexton
X = reshape(responses, size(responses,1)* size(responses,2), bankNum);
dis2textons = dist2(textons, X);
[~, indxtexton] = max(dis2textons);
feattexton = reshape(indxtexton, row, col);

% 图像边界处理
if winSize > 1
    colNumLeft = floor((winSize-1)/2);
    colNumRight = ceil((winSize-1)/2);
    for i=1:colNumLeft
        feattexton = [feattexton(:,1), feattexton];
        feattexton = [feattexton(1,:); feattexton];
    end
    for i=1:colNumRight
        feattexton = [feattexton, feattexton(:,size(feattexton, 2))];
        feattexton = [feattexton; feattexton(size(feattexton, 1), :)];
    end
else
    colNumLeft = 0;
    colNumRight = 0;
end

% 生成纹理柱状图 featIm
featIm = zeros(row, col, size(textons, 1));
for i=(1+colNumLeft):(size(feattexton, 1)-colNumRight)
    for j=(1+colNumLeft):(size(feattexton, 2)-colNumRight)
        window = feattexton((i-colNumLeft):(i+colNumRight), (j-colNumLeft):(j+colNumRight));
        frequency = tabulate(window(:));
        for k=1:size(frequency, 1)
            textonIndex = int64(frequency(k, 1));
            count = int32(frequency(k, 2));
            featIm(i-colNumLeft, j-colNumLeft, textonIndex) = featIm(i-colNumLeft, j-colNumLeft, textonIndex)+count;
        end
    end
end

return;

% 给定一张灰度图像，一个滤波器组，一份纹理编码集，构建纹理柱状图。
% 对于每个像素，基于每个纹理在其邻近范围（定义在固定大小的winSize 内的局
% 部窗口）内出现的频率。其中，texton是k*d 的矩阵。

计算两种图像分割结果

函数说明

给定一个原始图像为h*w*3 的RGB 彩色图像，计算两种图像分割的结果：一个是基于颜色特征的，另外一个是基于纹理特征。基于颜色特征的图像分割采用基于k-means 聚类的算法，其中用于聚类的颜色信息应该是出现在给定的图片中的。而基于纹理特征的图像分割应基于图像纹理基元直方图的k-means 聚类算法。其中：colorLabelIm 和textureLabelIm 是h*w 的矩阵，分别表示基于颜色和基于纹理的分割区域的标签。而numColorRegion 和numTextRegions 分别表示上述两种特征类型指定的理想分割数目，其他的参数如同以上所定义的。函数形式：

function [colorLabelIm, textureLabelIm] = compareSegmentations(origIm, bank, textons, winSize, numColorRegions, numTextureRegions)

实现思路

调用上述方法，计算原始图像 origIm 的颜色分类隶属度矩阵 colorLabelIm，以及纹理基元分类隶属度矩阵 textureLabelIm。

代码展示

function [colorLabelIm, textureLabelIm] = ...
    compareSegmentations(origIm, bank, textons, winSize, ...
    numColorRegions, numTextureRegions)

origIm = im2double(origIm);

% 获取颜色集，生成基于颜色分割的标签矩阵
colordata = reshape(origIm, size(origIm, 1)*size(origIm, 2), size(origIm, 3));
opts = statset('Display','final','MaxIter',1000);
[~, colorCenter] = kmeans(colordata, numColorRegions, 'Options', opts);
colorLabelIm = quantizeFeats(origIm, colorCenter);

% 获取纹理柱状图，计算纹理特征，生成基于纹理分割的标签矩阵
featIm = extractTextonHists(rgb2gray(origIm), bank, textons, winSize);
featImData = reshape(featIm, size(featIm, 1)*size(featIm, 2), size(featIm, 3));
[~, textureCenter] = kmeans(featImData, numTextureRegions, 'Options', opts);
textureLabelIm = quantizeFeats(featIm, textureCenter);

return