最近论文在用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. Eddins 11 % Digital Image Processing Using MATLAB, Prentice-Hall, 2004 12 % $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 end 19 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 bit 31 % 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 end 42 % 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'即miss classification 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 end 13 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 end 9 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);
