哈工大深研院数字图像处理第二次大作业：水果自动识别（2）HSV空间聚类及SIFT算法目标识别

老规矩，直接贴报告~

Programe list:

Programe was developed in the condition of Windows aswell as Linux server, programming language is Matlab (www.mathworks.com).

Classify.m, Kmeans.m: function for K-means clustering.

main_Kmeans.m: main function for K-means clustering.

HSVfeatureExtraction.m: extracting the HSV features for a picture.

HSVfeatureSorting.m: sorting the extracted HSV features.

HSVMED.m,HSVMED_accuracy.m: making the MEDclassifcation in hsv space and calculate its classfication accuracy.

sift.m,match.m, showkeys.m: function for SIFTalgorithm.

MainG1 G3 G5.m: making recognition among specific groupsand calculating its recognition accuracy (by 100 times examine).

Methodology:

The 100 images are from 10 kinds of fruits and each fruit contains 10images. First we dicided to use SIFT algorithm for all of the fruitsrecognition but the images are too big to process, even on super computerserver provided by Sun Yat-sen University. Hence,
we decided to make an advanceclassification by colors then apply SIFT among the same groups.

Full algorithm flow diagram is as follows:



Fig1. Algorithm flowdiagram
Firstly, perior to the color feature extraction, region of interest (ROI)extraction is the basic step for image processing, here we applied an edgedetection method called watershed. Algorithm of watershed is a segement methodbased on topology theory. The
basic concept is regarding the image as alandform of a map, altitude is defined by the gray level at that pixel. Local minimalvalue and its related region called catchment basin and its boundry calledwatershed. More specific, sketch map was shown in Fig.2:



Fig.2 Sketch map ofwatershed
To calculate thewatershed in an image, iterative algorithm should be applied; core algorithmsare sorting and submerging. Firstly gray levels of each pixel are sort in theforward direction. Then ‘wate’r was ‘poured’ on the image and watershed
willappear. ROI extraction was taken place in Hue of HSV space and the extractionresult is shown in Fig.3:



Fig.3 ROI extractionresult of watershed method
As we can see in Fig.3,watershed method can derive feasible results.
RGB color space was firstly applied for the color classfication feature,however, result by this feature was bad and then we decided to do thecalssfication in HSV space.

By extracting the HSV features, we can draw these features of eachpicture in HSV color space and calculated the mean feature of each group:



Fig.4 HSV features of imagesand groups in HSV space
As is shown in Fig.4, redpoints represent the HSV feature of each image and blue points represent theHSV feature of each group. We can conclude that pictures in same group are incorrelation with each other and have difference from pictures form
differentgroups. Then MED was applied for each images and try to assign them intospecific group.



Fig.5 Classifcation withMED in HSV color space
However, finalclassification result was not feasible, only 67% pictures are assigned in toright catalog. Then we decided to expand the size of groups and try to make theclassification perfect as well as less members in groups for a betterperformance
of post-SIFT process.
K-means clustering methodwas applied to those pictures, firstly the number of seed points was two and wecan get the correct result but each group has 40 and 60 pictures which cancause lots of trouble in further SIFT process. So the number of
seed pointsstarted to grow until clustering error happens, finally we can get five groupscontain different kinds of fruits:



Fig.6 Results of fruitsclassfication in HSV color space
After classification,Group2 and Group 4 contain only one kind of fruit which means classificationcan get 100% accuracy for the recognition of these two fruits, and Group1 andGroup3 contain three kinds of fruits and Group 5 contain 2 kinds of
fruits.
In the folder ‘/testsample’,G1-G4 only left one picture for giving the example and G5 are the ROIs ofdifferent two kinds of fruits which is ‘番石榴’ and ‘青苹果’, in order toindentify these two kinds of fruits of similar color, SIFT was applied.
SIFT is the algorithmwhich aims at detecting local features, by detecting the interest points of oneimage and descriptors of scale and orientation then making the match amongdifferent pictures, a feasible result can be obtained.



Fig.7 Flow diagram ofSIFT matching
By examing SIFT in thesame groups to obtain the accuracy, test samples were generated by random andother pictures are used for training sampls. After 100 time of random checkout,the highest identification probability can be 93%, which is acceptable
for ourdaily use.
Furter works can aim atthe better efficiency of ROI extraction and grouping, plastic bags are terriblefor texture feature extraction, testing on texture generated a bad result. Ifwe can remove the influence of plastic bags, I think texture features
will giveus some interesting results.
Reference:

《数字图像处理（Matlab版）》第二版 冈萨雷斯
Lowe SIFT 原文：
http://www.cs.ubc.ca/~lowe/papers/ijcv04.pdf
MATLAB 应用Sift算子的模式识别方法：/article/1820998.html
SIFT算法讲解：
/article/1820926.html
K-means 算法讲解：
/article/5017131.html
分水岭分割算法：
/article/8150158.html
贴代码：

%Name: HSV feature、 extraction
%Function: Extract ROI feature from HSV images
%Author:    Changle Zhang, zhangchanglehit@163.com
clc;
clear all;
close all;
%Initilization

pathname='E:\Pro2_15S158746_张常乐\Fruit Samples For Project2\雪山白苹果';           %修改文件夹名称
cd(pathname);

dirs=dir([,'*.jpg']);
dircell=struct2cell(dirs)';
filenames=dircell(:,1);
hsvmean=zeros(10,3);             %预分配平均参数
for number=1:10

rgbimg=imread(char(filenames(number)));
r=rgbimg(:,:,1);
g=rgbimg(:,:,2);
b=rgbimg(:,:,3);

grayimg=rgb2gray(rgbimg);

hsvimg=rgb2hsv(rgbimg);
H=hsvimg(:,:,1);
S=hsvimg(:,:,2);
V=hsvimg(:,:,3);

%rgb = imread('pears.png');%读取原图像
I = S;%转化为灰度图像
hy = fspecial('sobel');%sobel算子
hx = hy';
Iy = imfilter(double(I), hy, 'replicate');%滤波求y方向边缘
Ix = imfilter(double(I), hx, 'replicate');%滤波求x方向边缘
gradmag = sqrt(Ix.^2 + Iy.^2);%求摸
%3.分别对前景和背景进行标记：本例中使用形态学重建技术对前景对象进行标记，首先使用开操作，开操作之后可以去掉一些很小的目标。
se = strel('disk', 20);%圆形结构元素
Io = imopen(I, se);%形态学开操作
Ie = imerode(I, se);%对图像进行腐蚀
Iobr = imreconstruct(Ie, I);%形态学重建
Ioc = imclose(Io, se);%形态学关操作
Iobrd = imdilate(Iobr, se);%对图像进行膨胀
Iobrcbr = imreconstruct(imcomplement(Iobrd), imcomplement(Iobr));%形态学重建
Iobrcbr = imcomplement(Iobrcbr);%图像求反
bw = im2bw(Iobrcbr, graythresh(Iobrcbr));%转化为二值图像
%figure;
%imshow(bw), %显示二值图像
%title('Thresholded opening-closing by reconstruction')
bww=double(bw);
hsvimg(:,:,1)=bww.*hsvimg(:,:,1);
hsvimg(:,:,2)=bww.*hsvimg(:,:,2);
hsvimg(:,:,3)=bww.*hsvimg(:,:,3);
%figure;
%imshow(rgbimg);

imgsize=size(grayimg);
counter=0;
hsv=[0 0 0];
hsv=double(hsv);
for i=1:imgsize(1)
for j=1:imgsize(2)
if bww(i,j)>0
hsv(1)=hsv(1)+double(hsvimg(i,j,1));
hsv(2)=hsv(2)+double(hsvimg(i,j,2));
hsv(3)=hsv(3)+double(hsvimg(i,j,3));
counter=counter+1;
end
end
end
hsvmean(number,:)=hsv/counter;
end
allhsvmean=mean(hsvmean);
save('hsvmean.txt','hsvmean','-ascii');             %保存每幅图像的ROI内hsv均值参数信息
save('allhsvmean.txt','allhsvmean','-ascii');       %保存整个类内ROI hsv均值参数信息

%Name: HSV feature MED
%Function: Get the accury for a HSV MED results
%Author:    Changle Zhang, zhangchanglehit@163.com
clc;
clear all;
close all;
%Initilization
load chsvfeature.mat;

pathname='E:\cl\DIP\Pro2_15S158746_张常乐\Pro2_15S158746_张常乐\Fruit Samples For Project2';           %^^change this to make analysis among NC or PSD/or AD
cd(pathname);
dirs=dir;
dircell=struct2cell(dirs)';
filename=dircell(3:12,1);
results=cell(11,10);
for num=1:10
cd(char(filename(num)));
results(11,num)=filename(num);
hsvs=load('hsvmean.txt');
for i=1:10
test=hsvs(i,:);
distance=zeros(10,1);
for j=1:10
distance(j)=(test(1)-features(j,1))^2+(test(2)-features(j,2))^2+(test(3)-features(j,3))^2;
end
position=find(distance==min(distance));
results(i,num)=filenames(position);
end
cd(pathname);
end
counter=0;
for i=1:10
for j=1:10
if sum(size(char(results(1,i)))==size(char(results(j,i))))==2
if sum(size(char(results(1,i)))==size(char(results(j,i))))~=0
counter=counter+1;
end
end
end
end
accury=counter/100;
%save ('hsvresults.txt','results');

%Name: mainG5
%Function: Get the accury for SIFT in group 5
%Author:    Changle Zhang, zhangchanglehit@163.com
clc;
clear all;
close all;
%初始化
pathname='E:\cl\DIP\Pro2_15S158746_张常乐\Pro2_15S158746_张常乐\testsample\G5';           %^^change this to make analysis
cd(pathname);
dirs=dir([,'*.jpg']);
dircell=struct2cell(dirs)';
filenames=dircell(:,1);
accury=0;

for i=1:100
testsample=round(20*rand);

test=char(filenames(testsample));

alla=0;
allb=0;

for num=1:10
if num~=testsample
[numa(num),m1,m2]=match(test,char(filenames(num)));
alla=alla+m2;
end
end
a=sum(numa)/alla;

for num=11:20
if num~=testsample
[numb(num),m1,m2]=match(test,char(filenames(num)));
allb=allb+m2;
end
end

b=sum(numb)/allb;
if a>=b
test='番石榴';
testnum=1;
end

if a<b
test='青苹果';
testnum=2;
end
if ((testnum==1)&&(testsample<=10))||((testnum==1)&&(testsample>=10))
accury=accury+1;
end
end
accury

% num = match(image1, image2)
%
% This function reads two images, finds their SIFT features, and
%   displays lines connecting the matched keypoints.  A match is accepted
%   only if its distance is less than distRatio times the distance to the
%   second closest match.
% It returns the number of matches displayed.
%
% Example: match('scene.pgm','book.pgm');

function [num,m1,m2] = match(image1, image2)

% Find SIFT keypoints for each image
[im1, des1, loc1,m1] = sift(image1);
[im2, des2, loc2,m2] = sift(image2);

% For efficiency in Matlab, it is cheaper to compute dot products between
%  unit vectors rather than Euclidean distances.  Note that the ratio of
%  angles (acos of dot products of unit vectors) is a close approximation
%  to the ratio of Euclidean distances for small angles.
%
% distRatio: Only keep matches in which the ratio of vector angles from the
%   nearest to second nearest neighbor is less than distRatio.
distRatio = 0.62;

% For each descriptor in the first image, select its match to second image.
des2t = des2';                          % Precompute matrix transpose
for i = 1 : size(des1,1)
dotprods = des1(i,:) * des2t;        % Computes vector of dot products
[vals,indx] = sort(acos(dotprods));  % Take inverse cosine and sort results

% Check if nearest neighbor has angle less than distRatio times 2nd.
if (vals(1) < distRatio * vals(2))
match(i) = indx(1);
else
match(i) = 0;
end
end

% Create a new image showing the two images side by side.
% im3 = appendimages(im1,im2);

% Show a figure with lines joining the accepted matches.
%figure('Position', [100 100 size(im3,2) size(im3,1)]);
% colormap('gray');
% imagesc(im3);
% hold on;
% cols1 = size(im1,2);
% for i = 1: size(des1,1)
%   if (match(i) > 0)
%     line([loc1(i,2) loc2(match(i),2)+cols1], ...
%          [loc1(i,1) loc2(match(i),1)], 'Color', 'c');
%   end
% end
% hold off;
num = sum(match > 0);
fprintf('Found %d matches.\n', num);
display(m2)

% [image, descriptors, locs] = sift(imageFile)
%
% This function reads an image and returns its SIFT keypoints.
%   Input parameters:
%     imageFile: the file name for the image.
%
%   Returned:
%     image: the image array in double format
%     descriptors: a K-by-128 matrix, where each row gives an invariant
%         descriptor for one of the K keypoints.  The descriptor is a vector
%         of 128 values normalized to unit length.
%     locs: K-by-4 matrix, in which each row has the 4 values for a
%         keypoint location (row, column, scale, orientation).  The
%         orientation is in the range [-PI, PI] radians.
%
% Credits: Thanks for initial version of this program to D. Alvaro and
%          J.J. Guerrero, Universidad de Zaragoza (modified by D. Lowe)

function [image, descriptors, locs,m] = sift(imageFile)

% Load image
image = rgb2gray(imread(imageFile));
% imshow(image,[]);%%%%%%%%%%%%%%%%%%%%%
% If you have the Image Processing Toolbox, you can uncomment the following
%   lines to allow input of color images, which will be converted to grayscale.
% if isrgb(image)
%    image = rgb2gray(image);
% end

[rows, cols] = size(image);
% Convert into PGM imagefile, readable by "keypoints" executable
f = fopen('tmp.pgm', 'w');
if f == -1
error('Could not create file tmp.pgm.');
end
fprintf(f, 'P5\n%d\n%d\n255\n', cols, rows);
fwrite(f, image', 'uint8');
fclose(f);

% Call keypoints executable
if isunix
command = '!./sift ';
else
command = '!siftWin32 ';
end
command = [command ' <tmp.pgm >tmp.key'];
eval(command);

% Open tmp.key and check its header
g = fopen('tmp.key', 'r');
if g == -1
error('Could not open file tmp.key.');
end
[header, count] = fscanf(g, '%d %d', [1 2]);
if count ~= 2
error('Invalid keypoint file beginning.');
end
num = header(1);
len = header(2);
if len ~= 128
error('Keypoint descriptor length invalid (should be 128).');
end

% Creates the two output matrices (use known size for efficiency)
locs = double(zeros(num, 4));
descriptors = double(zeros(num, 128));

% Parse tmp.key
for i = 1:num
[vector, count] = fscanf(g, '%f %f %f %f', [1 4]); %row col scale ori
if count ~= 4
error('Invalid keypoint file format');
end
locs(i, :) = vector(1, :);

[descrip, count] = fscanf(g, '%d', [1 len]);
if (count ~= 128)
error('Invalid keypoint file value.');
end
% Normalize each input vector to unit length
descrip = descrip / sqrt(sum(descrip.^2));
descriptors(i, :) = descrip(1, :);
end
m=size(descriptors,1);
fclose(g);

% showkeys(image, locs)
%
% This function displays an image with SIFT keypoints overlayed.
%   Input parameters:
%     image: the file name for the image (grayscale)
%     locs: matrix in which each row gives a keypoint location (row,
%           column, scale, orientation)

function showkeys(image, locs)

disp('Drawing SIFT keypoints ...');

% Draw image with keypoints
figure('Position', [50 50 size(image,2) size(image,1)]);
colormap('gray');
imagesc(image);
hold on;
imsize = size(image);
for i = 1: size(locs,1)
% Draw an arrow, each line transformed according to keypoint parameters.
TransformLine(imsize, locs(i,:), 0.0, 0.0, 1.0, 0.0);
TransformLine(imsize, locs(i,:), 0.85, 0.1, 1.0, 0.0);
TransformLine(imsize, locs(i,:), 0.85, -0.1, 1.0, 0.0);
end
hold off;

% ------ Subroutine: TransformLine -------
% Draw the given line in the image, but first translate, rotate, and
% scale according to the keypoint parameters.
%
% Parameters:
%   Arrays:
%    imsize = [rows columns] of image
%    keypoint = [subpixel_row subpixel_column scale orientation]
%
%   Scalars:
%    x1, y1; begining of vector
%    x2, y2; ending of vector
function TransformLine(imsize, keypoint, x1, y1, x2, y2)

% The scaling of the unit length arrow is set to approximately the radius
%   of the region used to compute the keypoint descriptor.
len = 6 * keypoint(3);

% Rotate the keypoints by 'ori' = keypoint(4)
s = sin(keypoint(4));
c = cos(keypoint(4));

% Apply transform
r1 = keypoint(1) - len * (c * y1 + s * x1);
c1 = keypoint(2) + len * (- s * y1 + c * x1);
r2 = keypoint(1) - len * (c * y2 + s * x2);
c2 = keypoint(2) + len * (- s * y2 + c * x2);

line([c1 c2], [r1 r2], 'Color', 'c');

%画图用
load('hsvfeature.mat');
a=features;
plot3(a(:,1),a(:,2),a(:,3),'b.');
hold on

pathname='E:\cl\DIP\Pro2_15S158746_张常乐\Pro2_15S158746_张常乐\Fruit Samples For Project2';           %^^change this to make analysis among NC or PSD/or AD
cd(pathname);
dirs=dir;
dircell=struct2cell(dirs)';
filenames=dircell(3:12,1);
results=cell(11,10);
for num=1:10
cd(char(filenames(num)));
hsvs=load('hsvmean.txt');
plot3(hsvs(:,1),hsvs(:,2),hsvs(:,3),'r.');
hold on
cd(pathname);
end