www.gusucode.com > images 案例代码 matlab源码程序 > images/TextureSegmentationUsingGaborFiltersExample.m

    %% Texture Segmentation Using Gabor Filters
% This example shows how to use texture segmentation to identify regions
% based on their texture. The goal is to segment the dog from the bathroom
% floor. The segmentation is visually obvious because of the difference in
% texture between the regular, periodic pattern of the bathroom floor, and
% the regular, smooth texture of the dog's fur.
%
% From experimention, it is known that Gabor filters are a reasonable model
% of simple cells in the Mammalian vision system. Because of this, Gabor
% filters are throught to be a good model of how humans distinguish
% texture, and are therefore a useful model to use when designing
% algorithms to recognize texture. This example uses the basic approach
% described in (A. K. Jain and F. Farrokhnia, "Unsupervised Texture
% Segmentation Using Gabor Filters",1991) to perform texture segmentation.
%
% This example requires Statistics and Machine Learning Toolbox(TM).
% 
% Copyright 2015 The MathWorks, Inc.
%% Read and display input image
% Read and display the input image. This example shrinks the image to make
% example run more quickly.
A = imread('kobi.png');
A = imresize(A,0.25);
Agray = rgb2gray(A);
figure
imshow(A)
%% Design Array of Gabor Filters 
% Design an array of Gabor Filters which are tuned to different frequencies
% and orientations. The set of frequencies and orientations is designed to
% localize different, roughly orthogonal, subsets of frequency and
% orientation information in the input image. Regularly sample orientations
% between [0,150] degrees in steps of 30 degrees. Sample wavlength in
% increasing powers of two starting from 4/sqrt(2) up to the hypotenuse
% length of the input image. These combinations of frequency and
% orientation are taken from [Jain,1991] cited in the introduction.
imageSize = size(A);
numRows = imageSize(1);
numCols = imageSize(2);

wavelengthMin = 4/sqrt(2);
wavelengthMax = hypot(numRows,numCols);
n = floor(log2(wavelengthMax/wavelengthMin));
wavelength = 2.^(0:(n-2)) * wavelengthMin;

deltaTheta = 45;
orientation = 0:deltaTheta:(180-deltaTheta);

g = gabor(wavelength,orientation);
%% 
% Extract gabor magnitude features from source image. When working with
% Gabor filters, it is common to work with the magnitude response of each
% filter. Gabor magnitude response is also sometimes referred to as "Gabor
% Energy". Each MxN Gabor magnitude output image in gabormag(:,:,ind) is
% the output of the corresponding gabor filter g(ind).
gabormag = imgaborfilt(Agray,g);
%% Post-process the Gabor Magnitude Images into Gabor Features. 
% To use Gabor magnitude responses as features for use in classification,
% some post-processing is required. This post processing includes Gaussian
% smoothing, adding additional spatial information to the feature set,
% reshaping our feature set to the form expected by the |pca| and |kmeans|
% functions, and normalizing the feature information to a common variance
% and mean. 
%
% Each Gabor magnitude image contains some local variations, even within
% well segmented regions of constant texture. These local variations will
% throw off the segmentation. We can compensate for these variations using
% simple Gaussian low-pass filtering to smooth the Gabor magnitude
% information. We choose a sigma that is matched to the Gabor filter that
% extracted each feature. We introduce a smoothing term K that controls how
% much smoothing is applied to the Gabor magnitude responses.
for i = 1:length(g)
    sigma = 0.5*g(i).Wavelength;
    K = 3;
    gabormag(:,:,i) = imgaussfilt(gabormag(:,:,i),K*sigma); 
end
%%
% When constructing Gabor feature sets for classification, it is useful to
% add a map of spatial location information in both X and Y. This
% additional information allows the classifier to prefer groupings which
% are close together spatially.
X = 1:numCols;
Y = 1:numRows;
[X,Y] = meshgrid(X,Y);
featureSet = cat(3,gabormag,X);
featureSet = cat(3,featureSet,Y);
%%
% Reshape data into a matrix X of the form expected by the |kmeans|
% function. Each pixel in the image grid is a separate datapoint, and each
% plane in the variable |featureSet| is a separate feature. In this
% example, there is a separate feature for each filter in the Gabor filter
% bank, plus two additional features from the spatial information that was
% added in the previous step. In total, there are 24 Gabor features and 2
% spatial features for each pixel in our input image.
numPoints = numRows*numCols;
X = reshape(featureSet,numRows*numCols,[]);
%%
% Normalize features to be zero mean, unit variance.
X = bsxfun(@minus, X, mean(X));
X = bsxfun(@rdivide,X,std(X));

%% 
% Visualize feature set. To get a sense of what the Gabor magnitude
% features look like, Principal Component Analysis can be used to move from
% a 26-D representation of each pixel in the input image into a 1-D
% intensity value for each pixel.
coeff = pca(X);
feature2DImage = reshape(X*coeff(:,1),numRows,numCols);
figure
imshow(feature2DImage,[])
%%
% It is apparent in this visualization that there is sufficient variance in
% the Gabor feature information to obtain a good segmentation for this
% image. The dog is very dark compared to the floor because of the texture
% differences between the dog and the floor.
%% Classify Gabor Texture Features using kmeans
% Repeat k-means clustering five times to avoid local minima when searching
% for means that minimize objective function. The only prior information
% assumed in this example is how many distinct regions of texture are
% present in the image being segmented. There are two distinct regions in
% this case. This part of the example requires the Statistics and Machine
% Learning Toolbox.
L = kmeans(X,2,'Replicates',5);
%% 
% Visualize segmentation using |label2rgb|.
L = reshape(L,[numRows numCols]);
figure
imshow(label2rgb(L))
%% 
% Visualize the segmented image using |imshowpair|. Examine the foreground
% and background images that result from the mask BW that is associated
% with the label matrix L.
Aseg1 = zeros(size(A),'like',A);
Aseg2 = zeros(size(A),'like',A);
BW = L == 2;
BW = repmat(BW,[1 1 3]);
Aseg1(BW) = A(BW);
Aseg2(~BW) = A(~BW);
figure
imshowpair(Aseg1,Aseg2,'montage');