用MATLAB进行SVM分类

  最近论文在用SVM进行分类,目的是检测缺陷。缺陷有三种分别是孔洞,刮擦和划痕缺陷。 我用过libsvm和ddtools还有就是matlab中的svm函数 (svmtrain和svmclsassify),libsvm原来用的效果不好,我现在又忘了怎么用了,改天再把它捡起来吧,现在为应付毕业用的是 matlab的函数。 
       进行svm分类有这么几个步骤,第一是采集样本;提取样本的特征;进行cross validation验证支持向量机中的参数(一般都是选择的rbf的非线性支持向量机,参数有2个,即惩罚系数和rbf核的半径sigma)。
       1、在提取样本时,必须将样本图像缩放成统一大小的二值化(或灰度图像),因为我检测的缺陷中光照变化比较强烈,所以我只是选择了二值化图像,每个图像缩放成50*50的尺寸。注意不同大小的图像后面进行特征提取是没有意义的。
       2、其次提样本特征,对于缺陷图像分类,一般来说提取的都是特征矩,我提取的是图像的不变矩,它不随图像的旋转和缩放而产生变化(不变矩最后计算出来有7种参数),当然也可以提取多种矩然后用PCA的进行降维,matlab的PCA很容易,以前的博客上有介绍。
用matlab提取不变矩的算法网上有,叫做invmoments。
function phi = invmoments(F)
%INVMOMENTS Compute invariant moments of image.
%   PHI = INVMOMENTS(F) computes the moment invariants of the image
%   F. PHI is a seven-element row vector containing the moment
%   invariants as defined in equations (11.3-17) through (11.3-23) of
%   Gonzalez and Woods, Digital Image Processing, 2nd Ed.
%
%   F must be a 2-D, real, nonsparse, numeric or logical matrix.

%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.5 $  $Date: 2003/11/21 14:39:19 $

if (ndims(F) ~= 2) | issparse(F) | ~isreal(F) | ~(isnumeric(F) | ...
                                                    islogical(F))
   error(['F must be a 2-D, real, nonsparse, numeric or logical ' ...
          'matrix.']); 
end

F = double(F);

phi = compute_phi(compute_eta(compute_m(F)));
  
%-------------------------------------------------------------------%
function m = compute_m(F)

[M, N] = size(F);
[x, y] = meshgrid(1:N, 1:M);
  
% Turn x, y, and F into column vectors to make the summations a bit
% easier to compute in the following.
x = x(:);
y = y(:);
F = F(:);
  
% DIP equation (11.3-12)
m.m00 = sum(F);
% Protect against divide-by-zero warnings.
if (m.m00 == 0)
   m.m00 = eps;
end
% The other central moments:  
m.m10 = sum(x .* F);
m.m01 = sum(y .* F);
m.m11 = sum(x .* y .* F);
m.m20 = sum(x.^2 .* F);
m.m02 = sum(y.^2 .* F);
m.m30 = sum(x.^3 .* F);
m.m03 = sum(y.^3 .* F);
m.m12 = sum(x .* y.^2 .* F);
m.m21 = sum(x.^2 .* y .* F);

%-------------------------------------------------------------------%
function e = compute_eta(m)

% DIP equations (11.3-14) through (11.3-16).

xbar = m.m10 / m.m00;
ybar = m.m01 / m.m00;

e.eta11 = (m.m11 - ybar*m.m10) / m.m00^2;
e.eta20 = (m.m20 - xbar*m.m10) / m.m00^2;
e.eta02 = (m.m02 - ybar*m.m01) / m.m00^2;
e.eta30 = (m.m30 - 3 * xbar * m.m20 + 2 * xbar^2 * m.m10) / m.m00^2.5;
e.eta03 = (m.m03 - 3 * ybar * m.m02 + 2 * ybar^2 * m.m01) / m.m00^2.5;
e.eta21 = (m.m21 - 2 * xbar * m.m11 - ybar * m.m20 + ...
           2 * xbar^2 * m.m01) / m.m00^2.5;
e.eta12 = (m.m12 - 2 * ybar * m.m11 - xbar * m.m02 + ...
           2 * ybar^2 * m.m10) / m.m00^2.5;

%-------------------------------------------------------------------% 
function phi = compute_phi(e)

% DIP equations (11.3-17) through (11.3-23).

phi(1) = e.eta20 + e.eta02;
phi(2) = (e.eta20 - e.eta02)^2 + 4*e.eta11^2;
phi(3) = (e.eta30 - 3*e.eta12)^2 + (3*e.eta21 - e.eta03)^2;
phi(4) = (e.eta30 + e.eta12)^2 + (e.eta21 + e.eta03)^2;
phi(5) = (e.eta30 - 3*e.eta12) * (e.eta30 + e.eta12) * ...
         ( (e.eta30 + e.eta12)^2 - 3*(e.eta21 + e.eta03)^2 ) + ...
         (3*e.eta21 - e.eta03) * (e.eta21 + e.eta03) * ...
         ( 3*(e.eta30 + e.eta12)^2 - (e.eta21 + e.eta03)^2 );
phi(6) = (e.eta20 - e.eta02) * ( (e.eta30 + e.eta12)^2 - ...
                                 (e.eta21 + e.eta03)^2 ) + ...
         4 * e.eta11 * (e.eta30 + e.eta12) * (e.eta21 + e.eta03);
phi(7) = (3*e.eta21 - e.eta03) * (e.eta30 + e.eta12) * ...
         ( (e.eta30 + e.eta12)^2 - 3*(e.eta21 + e.eta03)^2 ) + ...
         (3*e.eta12 - e.eta30) * (e.eta21 + e.eta03) * ...
         ( 3*(e.eta30 + e.eta12)^2 - (e.eta21 + e.eta03)^2 );
phi(1)=abs(log(phi(1)));
phi(2)=abs(log(phi(2)));
phi(3)=abs(log(phi(3)));
phi(4)=abs(log(phi(4)));
phi(5)=abs(log(phi(5)));
phi(6)=abs(log(phi(6)));
phi(7)=abs(log(phi(7)));

提取纹理矩的的函数也有,叫做statxture.m  
       提取完7个不变矩后,将n个样本的不变矩组成n*7的矩阵,(此时矩阵中包括多个分类的样本)。注意,然后要对n*7的矩阵进行rescale,即每一列除以该列中的最大元素,变成[0,1]的值。
经过处理后假设样本为sample,标签为grp,此时就可以进行cross validation. 交叉验证使用的matlab函数为crossval , crossval可以用于分类验证和回归验证,其中分类验证的参数名为'mcr'即misclassification rate误分率,回归验证的参数名为 'mse'即minimum square error。crossval的函数原型为: 
             

mcr=crossval(‘mcr’,样本集,grp,’predfun’,@函数句柄,’partition’, cp)

                                                                                                                 

误分率,                                 标签                   分类函数模型      cvpartition类中对象,

要求最小值                                                                                    用于k-fold分支交叉验证

在对svm进行训练时,为提高训练模型的准确率,减少过学习可能性,常使用k-fold进行验证,所谓的k-fold就是将样本(这里指包含所有类的训练样本)分成k部分,轮流拿其中1部分作为test,拿其它n-1部分作为train,在matlab 中函数cvpartition负责进行分类,在我的程序中我将样本分为5类。
调用crossval就是要找到参数C(惩罚系数)和sigma,使误分率mcr最小。相应的matlab代码为:

%进行非线性分类以及使用cross-validation 

function yfit = ...

    crossfun(xtrain,ytrain,xtest,rbf_sigma,boxconstraint)

% Train the model on xtrain, ytrain, 

% and get predictions of class of xtest

svmStruct = svmtrain(xtrain,ytrain,'Kernel_Function','rbf',...

   'rbf_sigma',rbf_sigma,'boxconstraint',boxconstraint);

yfit = svmclassify(svmStruct,xtest);

%上面求的是predfun中的函数句柄指向函数。


%下面求得是mcr误分率最小值 

%sample是样本矩,grp是标签

c=cvpartition(grp,’k’,5);  %我用5 分支

minfn = @(z)crossval('mcr',sample,grp,'Predfun', ...  %sample和grp换上你自己的样本和标签

    @(xtrain,ytrain,xtest)crossfun(xtrain,ytrain,...

    xtest,exp(z(1)),exp(z(2))),'partition',c);

opts = optimset('TolX',5e-4,'TolFun',5e-4); TolX x parameter tolerance  function %value tolerance

%最后要找最小值  这行代码是示例,真正用fminsearch找最小值需要用网格尽心定位的

[searchmin fval] = fminsearch(minfn,randn(2,1),opts) 这里是用fminsearch找最小值fval

在台湾大学林智仁那边,他提倡采用网格法找到最佳参数,我按他的方法定义了粗细两重网格:

C=[-5 -3 -1 1 3 5 7 9 11 13 15]; 

rbf_sigma=[-15 -13 -11 -9 -7 -5 -3 -1 1 3];

c=cvpartition(grp,'k',5); 

minfn = @(z)crossval('mcr',sample,grp,'Predfun'...

    @(xtrain,ytrain,xtest)crossfun(xtrain,ytrain,...

    xtest,exp(z(1)),exp(z(2))),'partition',c);

opts = optimset('TolX',5e-4,'TolFun',5e-4); 

for i=1:size(C,2)

    for j=1:size(rbf_sigma,2)

        [searchmin fval] = fminsearch(minfn,[rbf_sigma(j);C(i)],opts);

k(i,j)=fval;

    end

end

然后再通过找到k(i,j)最小值时的C和rbf_sigma的邻近区域组成细网格

C1=5:0.5:11;

sigma1=-1:0.25:3; 

k1(1:size(C1,2),1:size(sigma1,2))=zeros(size(C1,2),size(sigma1,2));

for i=1:size(C1,2)

    for j=1:size(sigma1,2)

        [searchmin fval] = fminsearch(minfn,[sigma1(1,j);C1(1,i)],opts);

        k1(i,j)=fval;

    end

end

当确定了C和sigma的值后,将sample,grp,C,rbf_sigma的值代入matlab程序svmtrain,求出svmStruct结构体,里面包含了支持向量等参数

svmStruct= svmtrain(xtrain,ytrain,'Kernel_Function','rbf',...
   'rbf_sigma',sigma,'boxconstraint',C);

把svmStruct存起来,后面用的时候load svmStruct就行了。 


对新来的缺陷图像进行分类时,首先把它缩成50*50的二值化图像,然后调用invmoments提取7个不变矩,然后要除以原来训练样本中得到的每一列元素的最大值进行归一化(记得前面我们说了要对不变矩7列中的每一列进行rescale吧,同时每一列的最大值还得保留起来,现在还要用)

然后就方便了

yfit = svmclassify(svmStruct, test);

 就可得到分类结果了。

  • 4
    点赞
  • 35
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 3
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 3
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

怀想天空2010

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值