首页 > 代码库 > 缺陷检测之SVM篇（update）
缺陷检测之SVM篇（update）
最近论文在用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。
1 function phi = invmoments(F) 2 %INVMOMENTS Compute invariant moments of image. 3 % PHI = INVMOMENTS(F) computes the moment invariants of the image 4 % F. PHI is a seven-element row vector containing the moment 5 % invariants as defined in equations (11.3-17) through (11.3-23) of 6 % Gonzalez and Woods, Digital Image Processing, 2nd Ed. 7 % 8 % F must be a 2-D, real, nonsparse, numeric or logical matrix. 9 10 % Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins11 % Digital Image Processing Using MATLAB, Prentice-Hall, 200412 % $Revision: 1.5 $ $Date: 2003/11/21 14:39:19 $13 14 if (ndims(F) ~= 2) | issparse(F) | ~isreal(F) | ~(isnumeric(F) | ...15 islogical(F))16 error([‘F must be a 2-D, real, nonsparse, numeric or logical ‘ ...17 ‘matrix.‘]); 18 end19 20 F = double(F);21 22 phi = compute_phi(compute_eta(compute_m(F)));23 24 %-------------------------------------------------------------------%25 function m = compute_m(F)26 27 [M, N] = size(F);28 [x, y] = meshgrid(1:N, 1:M);29 30 % Turn x, y, and F into column vectors to make the summations a bit31 % easier to compute in the following.32 x = x(:);33 y = y(:);34 F = F(:);35 36 % DIP equation (11.3-12)37 m.m00 = sum(F);38 % Protect against divide-by-zero warnings.39 if (m.m00 == 0)40 m.m00 = eps;41 end42 % The other central moments: 43 m.m10 = sum(x .* F);44 m.m01 = sum(y .* F);45 m.m11 = sum(x .* y .* F);46 m.m20 = sum(x.^2 .* F);47 m.m02 = sum(y.^2 .* F);48 m.m30 = sum(x.^3 .* F);49 m.m03 = sum(y.^3 .* F);50 m.m12 = sum(x .* y.^2 .* F);51 m.m21 = sum(x.^2 .* y .* F);52 53 %-------------------------------------------------------------------%54 function e = compute_eta(m)55 56 % DIP equations (11.3-14) through (11.3-16).57 58 xbar = m.m10 / m.m00;59 ybar = m.m01 / m.m00;60 61 e.eta11 = (m.m11 - ybar*m.m10) / m.m00^2;62 e.eta20 = (m.m20 - xbar*m.m10) / m.m00^2;63 e.eta02 = (m.m02 - ybar*m.m01) / m.m00^2;64 e.eta30 = (m.m30 - 3 * xbar * m.m20 + 2 * xbar^2 * m.m10) / m.m00^2.5;65 e.eta03 = (m.m03 - 3 * ybar * m.m02 + 2 * ybar^2 * m.m01) / m.m00^2.5;66 e.eta21 = (m.m21 - 2 * xbar * m.m11 - ybar * m.m20 + ...67 2 * xbar^2 * m.m01) / m.m00^2.5;68 e.eta12 = (m.m12 - 2 * ybar * m.m11 - xbar * m.m02 + ...69 2 * ybar^2 * m.m10) / m.m00^2.5;70 71 %-------------------------------------------------------------------% 72 function phi = compute_phi(e)73 74 % DIP equations (11.3-17) through (11.3-23).75 76 phi(1) = e.eta20 + e.eta02;77 phi(2) = (e.eta20 - e.eta02)^2 + 4*e.eta11^2;78 phi(3) = (e.eta30 - 3*e.eta12)^2 + (3*e.eta21 - e.eta03)^2;79 phi(4) = (e.eta30 + e.eta12)^2 + (e.eta21 + e.eta03)^2;80 phi(5) = (e.eta30 - 3*e.eta12) * (e.eta30 + e.eta12) * ...81 ( (e.eta30 + e.eta12)^2 - 3*(e.eta21 + e.eta03)^2 ) + ...82 (3*e.eta21 - e.eta03) * (e.eta21 + e.eta03) * ...83 ( 3*(e.eta30 + e.eta12)^2 - (e.eta21 + e.eta03)^2 );84 phi(6) = (e.eta20 - e.eta02) * ( (e.eta30 + e.eta12)^2 - ...85 (e.eta21 + e.eta03)^2 ) + ...86 4 * e.eta11 * (e.eta30 + e.eta12) * (e.eta21 + e.eta03);87 phi(7) = (3*e.eta21 - e.eta03) * (e.eta30 + e.eta12) * ...88 ( (e.eta30 + e.eta12)^2 - 3*(e.eta21 + e.eta03)^2 ) + ...89 (3*e.eta12 - e.eta30) * (e.eta21 + e.eta03) * ...90 ( 3*(e.eta30 + e.eta12)^2 - (e.eta21 + e.eta03)^2 );91 phi(1)=abs(log(phi(1)));92 phi(2)=abs(log(phi(2)));93 phi(3)=abs(log(phi(3)));94 phi(4)=abs(log(phi(4)));95 phi(5)=abs(log(phi(5)));96 phi(6)=abs(log(phi(6)));97 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代码为：
1 %进行非线性分类以及使用cross-validation 2 function yfit = ... 3 crossfun(xtrain,ytrain,xtest,rbf_sigma,boxconstraint) 4 % Train the model on xtrain, ytrain, 5 % and get predictions of class of xtest 6 svmStruct = svmtrain(xtrain,ytrain,‘Kernel_Function‘,‘rbf‘,... 7 ‘rbf_sigma‘,rbf_sigma,‘boxconstraint‘,boxconstraint); 8 yfit = svmclassify(svmStruct,xtest); 9 %上面求的是’predfun’中的函数句柄指向函数。10 11 %下面求得是mcr误分率最小值 12 %sample是样本矩，grp是标签13 c=cvpartition(grp,’k’,5); %我用5个 分支14 minfn = @(z)crossval(‘mcr‘,sample,grp,‘Predfun‘, ... %sample和grp换上你自己的样本和标签15 @(xtrain,ytrain,xtest)crossfun(xtrain,ytrain,...16 xtest,exp(z(1)),exp(z(2))),‘partition‘,c);17 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；在台湾大学林智仁那边，他提倡采用网格法找到最佳参数，我按他的方法定义了粗细两重网格：
1 C=[-5 -3 -1 1 3 5 7 9 11 13 15]; 2 rbf_sigma=[-15 -13 -11 -9 -7 -5 -3 -1 1 3]; 3 c=cvpartition(grp,‘k‘,5); 4 minfn = @(z)crossval(‘mcr‘,sample,grp,‘Predfun‘, ... 5 @(xtrain,ytrain,xtest)crossfun(xtrain,ytrain,... 6 xtest,exp(z(1)),exp(z(2))),‘partition‘,c); 7 opts = optimset(‘TolX‘,5e-4,‘TolFun‘,5e-4); 8 for i=1:size(C,2) 9 for j=1:size(rbf_sigma,2)10 [searchmin fval] = fminsearch(minfn,[rbf_sigma(j);C(i)],opts);11 k(i,j)=fval;12 end13 end
然后再通过找到k(i,j)最小值时的C和rbf_sigma的邻近区域组成细网格
1 C1=5:0.5:11;2 sigma1=-1:0.25:3; 3 k1(1:size(C1,2),1:size(sigma1,2))=zeros(size(C1,2),size(sigma1,2));4 for i=1:size(C1,2)5 for j=1:size(sigma1,2)6 [searchmin fval] = fminsearch(minfn,[sigma1(1,j);C1(1,i)],opts);7 k1(i,j)=fval;8 end9 end
当确定了C和sigma的值后，将sample,grp,C,rbf_sigma的值代入matlab程序svmtrain，求出svmStruct结构体，里面包含了支持向量等参数
1 svmStruct=svmtrain(xtrain,ytrain,‘Kernel_Function‘,‘rbf‘,...‘rbf_sigma‘,sigma,‘boxconstraint‘,C);
把svmStruct存起来，后面用的时候load svmStruct就行了。
对新来的缺陷图像进行分类时，首先把它缩成50*50的二值化图像，然后调用invmoments提取7个不变矩，然后要除以原来训练样本中得到的每一列元素的最大值进行归一化(记得前面我们说了要对不变矩7列中的每一列进行rescale吧，同时每一列的最大值还得保留起来，现在还要用)然后就方便了,就可得到分类结果了。
1 yfit = svmclassify(svmStruct, test);
缺陷检测之SVM篇（update）