26 #include "base/Optionpk.h" 27 #include "algorithms/ConfusionMatrix.h" 28 #include "algorithms/FeatureSelector.h" 30 #include "algorithms/CostFactorySVM.h" 31 #include "algorithms/svm.h" 32 #include "imageclasses/ImgReaderOgr.h" 105 #define Malloc(type,n) (type *)malloc((n)*sizeof(type)) 107 double objFunction(
const std::vector<double> &x, std::vector<double> &grad,
void *my_func_data);
110 map<string,short> classValueMap;
111 vector<std::string> nameVector;
112 vector<unsigned int> nctraining;
113 vector<unsigned int> nctest;
114 Optionpk<std::string> svm_type_opt(
"svmt",
"svmtype",
"type of SVM (C_SVC, nu_SVC,one_class, epsilon_SVR, nu_SVR)",
"C_SVC");
115 Optionpk<std::string> kernel_type_opt(
"kt",
"kerneltype",
"type of kernel function (linear,polynomial,radial,sigmoid) ",
"radial");
117 Optionpk<float> coef0_opt(
"c0",
"coef0",
"coef0 in kernel function",0);
118 Optionpk<float> nu_opt(
"nu",
"nu",
"the parameter nu of nu-SVC, one-class SVM, and nu-SVR",0.5);
119 Optionpk<float> epsilon_loss_opt(
"eloss",
"eloss",
"the epsilon in loss function of epsilon-SVR",0.1);
120 Optionpk<int> cache_opt(
"cache",
"cache",
"cache memory size in MB",100);
121 Optionpk<float> epsilon_tol_opt(
"etol",
"etol",
"the tolerance of termination criterion",0.001);
122 Optionpk<bool> shrinking_opt(
"shrink",
"shrink",
"whether to use the shrinking heuristics",
false);
123 Optionpk<bool> prob_est_opt(
"pe",
"probest",
"whether to train a SVC or SVR model for probability estimates",
true,2);
124 Optionpk<bool> costfunction_opt(
"cf",
"cf",
"use Overall Accuracy instead of kappa",
false);
128 Optionpk<short> classvalue_opt(
"r",
"reclass",
"list of class values (use same order as in class opt).");
129 Optionpk<short> verbose_opt(
"v",
"verbose",
"use 1 to output intermediate results for plotting",0,2);
131 double objFunction(
const std::vector<double> &x, std::vector<double> &grad,
void *my_func_data){
133 assert(grad.empty());
134 vector<Vector2d<float> > *tf=
reinterpret_cast<vector<Vector2d<float>
>*> (my_func_data);
137 double error=1.0/epsilon_tol_opt[0];
141 CostFactorySVM costfactory(svm_type_opt[0], kernel_type_opt[0], kernel_degree_opt[0], gamma, coef0_opt[0], ccost, nu_opt[0], epsilon_loss_opt[0], cache_opt[0], epsilon_tol_opt[0], shrinking_opt[0], prob_est_opt[0], cv_opt[0], verbose_opt[0]);
147 costfactory.setCv(cv_opt[0]);
149 if(classname_opt.size()){
150 assert(classname_opt.size()==classvalue_opt.size());
151 for(
int iclass=0;iclass<classname_opt.size();++iclass)
152 costfactory.setClassValueMap(classname_opt[iclass],classvalue_opt[iclass]);
155 costfactory.setNameVector(nameVector);
157 for(
int iname=0;iname<nameVector.size();++iname){
158 if(costfactory.getClassValueMap().empty()){
159 costfactory.pushBackClassName(nameVector[iname]);
162 else if(costfactory.getClassIndex(type2string<short>((costfactory.getClassValueMap())[nameVector[iname]]))<0)
163 costfactory.pushBackClassName(type2string<short>((costfactory.getClassValueMap())[nameVector[iname]]));
166 costfactory.setNcTraining(nctraining);
167 costfactory.setNcTest(nctest);
169 kappa=costfactory.getCost(*tf);
173 int main(
int argc,
char *argv[])
175 map<short,int> reclassMap;
176 vector<int> vreclass;
177 Optionpk<string> training_opt(
"t",
"training",
"training vector file. A single vector file contains all training features (must be set as: b0, b1, b2,...) for all classes (class numbers identified by label option).");
178 Optionpk<float> ccost_opt(
"cc",
"ccost",
"min and max boundaries the parameter C of C-SVC, epsilon-SVR, and nu-SVR (optional: initial value)",1);
179 Optionpk<float> gamma_opt(
"g",
"gamma",
"min max boundaries for gamma in kernel function (optional: initial value)",0);
180 Optionpk<double> stepcc_opt(
"stepcc",
"stepcc",
"multiplicative step for ccost in GRID search",2);
181 Optionpk<double> stepg_opt(
"stepg",
"stepg",
"multiplicative step for gamma in GRID search",2);
184 Optionpk<string> label_opt(
"label",
"label",
"identifier for class label in training vector file.",
"label");
186 Optionpk<unsigned int> balance_opt(
"bal",
"balance",
"balance the input data to this number of samples for each class", 0);
187 Optionpk<bool> random_opt(
"random",
"random",
"in case of balance, randomize input data",
true);
188 Optionpk<int> minSize_opt(
"min",
"min",
"if number of training pixels is less then min, do not take this class into account", 0);
189 Optionpk<unsigned short> band_opt(
"b",
"band",
"band index (starting from 0, either use band option or use start to end)");
192 Optionpk<double> offset_opt(
"offset",
"offset",
"offset value for each spectral band input features: refl[band]=(DN[band]-offset[band])/scale[band]", 0.0);
193 Optionpk<double> scale_opt(
"scale",
"scale",
"scale value for each spectral band input features: refl=(DN[band]-offset[band])/scale[band] (use 0 if scale min and max in each band to -1.0 and 1.0)", 0.0);
196 Optionpk<double> tolerance_opt(
"tol",
"tolerance",
"relative tolerance for stopping criterion",0.0001);
198 input_opt.setHide(1);
199 tlayer_opt.setHide(1);
200 label_opt.setHide(1);
201 balance_opt.setHide(1);
202 random_opt.setHide(1);
203 minSize_opt.setHide(1);
205 bstart_opt.setHide(1);
207 offset_opt.setHide(1);
208 scale_opt.setHide(1);
209 svm_type_opt.setHide(1);
210 kernel_type_opt.setHide(1);
211 kernel_degree_opt.setHide(1);
212 coef0_opt.setHide(1);
214 epsilon_loss_opt.setHide(1);
215 cache_opt.setHide(1);
216 epsilon_tol_opt.setHide(1);
217 shrinking_opt.setHide(1);
218 prob_est_opt.setHide(1);
220 costfunction_opt.setHide(1);
221 maxit_opt.setHide(1);
222 tolerance_opt.setHide(1);
224 classname_opt.setHide(1);
225 classvalue_opt.setHide(1);
229 doProcess=training_opt.retrieveOption(argc,argv);
230 ccost_opt.retrieveOption(argc,argv);
231 gamma_opt.retrieveOption(argc,argv);
232 stepcc_opt.retrieveOption(argc,argv);
233 stepg_opt.retrieveOption(argc,argv);
234 input_opt.retrieveOption(argc,argv);
235 tlayer_opt.retrieveOption(argc,argv);
236 label_opt.retrieveOption(argc,argv);
237 balance_opt.retrieveOption(argc,argv);
238 random_opt.retrieveOption(argc,argv);
239 minSize_opt.retrieveOption(argc,argv);
240 band_opt.retrieveOption(argc,argv);
241 bstart_opt.retrieveOption(argc,argv);
242 bend_opt.retrieveOption(argc,argv);
243 offset_opt.retrieveOption(argc,argv);
244 scale_opt.retrieveOption(argc,argv);
245 svm_type_opt.retrieveOption(argc,argv);
246 kernel_type_opt.retrieveOption(argc,argv);
247 kernel_degree_opt.retrieveOption(argc,argv);
248 coef0_opt.retrieveOption(argc,argv);
249 nu_opt.retrieveOption(argc,argv);
250 epsilon_loss_opt.retrieveOption(argc,argv);
251 cache_opt.retrieveOption(argc,argv);
252 epsilon_tol_opt.retrieveOption(argc,argv);
253 shrinking_opt.retrieveOption(argc,argv);
254 prob_est_opt.retrieveOption(argc,argv);
255 cv_opt.retrieveOption(argc,argv);
256 costfunction_opt.retrieveOption(argc,argv);
257 maxit_opt.retrieveOption(argc,argv);
258 tolerance_opt.retrieveOption(argc,argv);
260 classname_opt.retrieveOption(argc,argv);
261 classvalue_opt.retrieveOption(argc,argv);
262 verbose_opt.retrieveOption(argc,argv);
264 catch(
string predefinedString){
265 std::cout << predefinedString << std::endl;
270 cout <<
"Usage: pkoptsvm -t training" << endl;
272 std::cout <<
"short option -h shows basic options only, use long option --help to show all options" << std::endl;
276 assert(training_opt.size());
280 if(verbose_opt[0]>=1){
282 std::cout <<
"input filename: " << input_opt[0] << std::endl;
283 std::cout <<
"training vector file: " << std::endl;
284 for(
int ifile=0;ifile<training_opt.size();++ifile)
285 std::cout << training_opt[ifile] << std::endl;
286 std::cout <<
"verbose: " << verbose_opt[0] << std::endl;
289 unsigned int totalSamples=0;
290 unsigned int totalTestSamples=0;
292 unsigned short nclass=0;
296 vector<double> offset;
297 vector<double> scale;
298 vector< Vector2d<float> > trainingPixels;
299 vector< Vector2d<float> > testPixels;
315 if(bstart_opt.size()){
316 if(bend_opt.size()!=bstart_opt.size()){
317 string errorstring=
"Error: options for start and end band indexes must be provided as pairs, missing end band";
321 for(
int ipair=0;ipair<bstart_opt.size();++ipair){
322 if(bend_opt[ipair]<=bstart_opt[ipair]){
323 string errorstring=
"Error: index for end band must be smaller then start band";
326 for(
int iband=bstart_opt[ipair];iband<=bend_opt[ipair];++iband)
327 band_opt.push_back(iband);
332 cerr << error << std::endl;
337 std::sort(band_opt.begin(),band_opt.end());
340 if(classname_opt.size()){
341 assert(classname_opt.size()==classvalue_opt.size());
342 for(
int iclass=0;iclass<classname_opt.size();++iclass)
343 classValueMap[classname_opt[iclass]]=classvalue_opt[iclass];
348 vector<string> fields;
350 trainingPixels.clear();
352 map<string,Vector2d<float> > trainingMap;
353 map<string,Vector2d<float> > testMap;
354 if(verbose_opt[0]>=1)
355 std::cout <<
"reading training file " << training_opt[0] << std::endl;
359 totalSamples=trainingReader.readDataImageOgr(trainingMap,fields,band_opt,label_opt[0],tlayer_opt,verbose_opt[0]);
360 if(input_opt.size()){
362 totalTestSamples=inputReader.readDataImageOgr(testMap,fields,band_opt,label_opt[0],tlayer_opt,verbose_opt[0]);
367 totalSamples=trainingReader.readDataImageOgr(trainingMap,fields,0,0,label_opt[0],tlayer_opt,verbose_opt[0]);
368 if(input_opt.size()){
370 totalTestSamples=inputReader.readDataImageOgr(testMap,fields,0,0,label_opt[0],tlayer_opt,verbose_opt[0]);
373 trainingReader.close();
375 if(trainingMap.size()<2){
379 string errorstring=
"Error: could not read at least two classes from training input file";
382 if(input_opt.size()&&testMap.size()<2){
383 string errorstring=
"Error: could not read at least two classes from test input file";
388 cerr << error << std::endl;
392 cerr <<
"error caught" << std::endl;
402 std::cout <<
"training pixels: " << std::endl;
403 map<string,Vector2d<float> >::iterator mapit;
404 mapit=trainingMap.begin();
405 while(mapit!=trainingMap.end()){
406 if(classValueMap.size()){
408 if(classValueMap[mapit->first]>0){
410 std::cout << mapit->first <<
" -> " << classValueMap[mapit->first] << std::endl;
413 std::cerr <<
"Error: names in classname option are not complete, please check names in training vector and make sure classvalue is > 0" << std::endl;
418 if((mapit->second).size()<minSize_opt[0]){
419 trainingMap.erase(mapit);
422 nameVector.push_back(mapit->first);
423 trainingPixels.push_back(mapit->second);
425 std::cout << mapit->first <<
": " << (mapit->second).size() <<
" samples" << std::endl;
430 nclass=trainingPixels.size();
431 if(classname_opt.size())
432 assert(nclass==classname_opt.size());
433 nband=trainingPixels[0][0].size()-2;
435 mapit=testMap.begin();
436 while(mapit!=testMap.end()){
437 if(classValueMap.size()){
439 if(classValueMap[mapit->first]>0){
443 std::cerr <<
"Error: names in classname option are not complete, please check names in test vector and make sure classvalue is > 0" << std::endl;
448 testPixels.push_back(mapit->second);
450 std::cout << mapit->first <<
": " << (mapit->second).size() <<
" samples" << std::endl;
453 if(input_opt.size()){
454 assert(nclass==testPixels.size());
455 assert(nband=testPixels[0][0].size()-2);
461 if(balance_opt[0]>0){
465 for(
int iclass=0;iclass<nclass;++iclass){
466 if(trainingPixels[iclass].size()>balance_opt[0]){
467 while(trainingPixels[iclass].size()>balance_opt[0]){
468 int index=rand()%trainingPixels[iclass].size();
469 trainingPixels[iclass].erase(trainingPixels[iclass].begin()+index);
473 int oldsize=trainingPixels[iclass].size();
474 for(
int isample=trainingPixels[iclass].size();isample<balance_opt[0];++isample){
475 int index = rand()%oldsize;
476 trainingPixels[iclass].push_back(trainingPixels[iclass][index]);
479 totalSamples+=trainingPixels[iclass].size();
481 assert(totalSamples==nclass*balance_opt[0]);
486 offset.resize(nband);
488 if(offset_opt.size()>1)
489 assert(offset_opt.size()==nband);
490 if(scale_opt.size()>1)
491 assert(scale_opt.size()==nband);
492 for(
int iband=0;iband<nband;++iband){
494 std::cout <<
"scaling for band" << iband << std::endl;
495 offset[iband]=(offset_opt.size()==1)?offset_opt[0]:offset_opt[iband];
496 scale[iband]=(scale_opt.size()==1)?scale_opt[0]:scale_opt[iband];
499 float theMin=trainingPixels[0][0][iband+startBand];
500 float theMax=trainingPixels[0][0][iband+startBand];
501 for(
int iclass=0;iclass<nclass;++iclass){
502 for(
int isample=0;isample<trainingPixels[iclass].size();++isample){
503 if(theMin>trainingPixels[iclass][isample][iband+startBand])
504 theMin=trainingPixels[iclass][isample][iband+startBand];
505 if(theMax<trainingPixels[iclass][isample][iband+startBand])
506 theMax=trainingPixels[iclass][isample][iband+startBand];
509 offset[iband]=theMin+(theMax-theMin)/2.0;
510 scale[iband]=(theMax-theMin)/2.0;
511 if(verbose_opt[0]>1){
512 std::cout <<
"Extreme image values for band " << iband <<
": [" << theMin <<
"," << theMax <<
"]" << std::endl;
513 std::cout <<
"Using offset, scale: " << offset[iband] <<
", " << scale[iband] << std::endl;
514 std::cout <<
"scaled values for band " << iband <<
": [" << (theMin-offset[iband])/scale[iband] <<
"," << (theMax-offset[iband])/scale[iband] <<
"]" << std::endl;
526 if(verbose_opt[0]>=1){
527 std::cout <<
"number of bands: " << nband << std::endl;
528 std::cout <<
"number of classes: " << nclass << std::endl;
536 nctraining.resize(nclass);
537 nctest.resize(nclass);
538 vector< Vector2d<float> > trainingFeatures(nclass);
539 for(
int iclass=0;iclass<nclass;++iclass){
540 if(verbose_opt[0]>=1)
541 std::cout <<
"calculating features for class " << iclass << std::endl;
542 nctraining[iclass]=trainingPixels[iclass].size();
543 if(verbose_opt[0]>=1)
544 std::cout <<
"nctraining[" << iclass <<
"]: " << nctraining[iclass] << std::endl;
545 if(testPixels.size()>iclass){
546 nctest[iclass]=testPixels[iclass].size();
547 if(verbose_opt[0]>=1){
548 std::cout <<
"nctest[" << iclass <<
"]: " << nctest[iclass] << std::endl;
554 trainingFeatures[iclass].resize(nctraining[iclass]+nctest[iclass]);
555 for(
int isample=0;isample<nctraining[iclass];++isample){
557 for(
int iband=0;iband<nband;++iband){
558 assert(trainingPixels[iclass].size()>isample);
559 assert(trainingPixels[iclass][isample].size()>iband+startBand);
560 assert(offset.size()>iband);
561 assert(scale.size()>iband);
562 float value=trainingPixels[iclass][isample][iband+startBand];
563 trainingFeatures[iclass][isample].push_back((value-offset[iband])/scale[iband]);
567 for(
int isample=0;isample<nctest[iclass];++isample){
569 for(
int iband=0;iband<nband;++iband){
570 assert(testPixels[iclass].size()>isample);
571 assert(testPixels[iclass][isample].size()>iband+startBand);
572 assert(offset.size()>iband);
573 assert(scale.size()>iband);
574 float value=testPixels[iclass][isample][iband+startBand];
576 trainingFeatures[iclass][nctraining[iclass]+isample].push_back((value-offset[iband])/scale[iband]);
579 assert(trainingFeatures[iclass].size()==nctraining[iclass]+nctest[iclass]);
582 assert(ccost_opt.size()>1);
583 if(ccost_opt.size()<3)
584 ccost_opt.push_back(sqrt(ccost_opt[0]*ccost_opt[1]));
585 assert(gamma_opt.size()>1);
586 if(gamma_opt.size()<3)
587 gamma_opt.push_back(sqrt(gamma_opt[0]*gamma_opt[1]));
588 assert(ccost_opt.size()==3);
589 assert(gamma_opt.size()==3);
590 assert(gamma_opt[0]<gamma_opt[1]);
591 assert(gamma_opt[0]<gamma_opt[2]);
592 assert(gamma_opt[2]<gamma_opt[1]);
593 assert(ccost_opt[0]<ccost_opt[1]);
594 assert(ccost_opt[0]<ccost_opt[2]);
595 assert(ccost_opt[2]<ccost_opt[1]);
597 std::vector<double> x(2);
606 const char* pszMessage;
607 void* pProgressArg=NULL;
608 GDALProgressFunc pfnProgress=GDALTermProgress;
611 pfnProgress(progress,pszMessage,pProgressArg);
612 double ncost=log(ccost_opt[1])/log(stepcc_opt[0])-log(ccost_opt[0])/log(stepcc_opt[0]);
613 double ngamma=log(gamma_opt[1])/log(stepg_opt[0])-log(gamma_opt[0])/log(stepg_opt[0]);
614 for(
double ccost=ccost_opt[0];ccost<=ccost_opt[1];ccost*=stepcc_opt[0]){
615 for(
double gamma=gamma_opt[0];gamma<=gamma_opt[1];gamma*=stepg_opt[0]){
618 std::vector<double> theGrad;
620 kappa=objFunction(x,theGrad,&trainingFeatures);
627 std::cout << ccost <<
" " << gamma <<
" " << kappa<< std::endl;
628 progress+=1.0/ncost/ngamma;
630 pfnProgress(progress,pszMessage,pProgressArg);
635 pfnProgress(progress,pszMessage,pProgressArg);
686 std::cout <<
" --ccost " << x[0];
687 std::cout <<
" --gamma " << x[1];
688 std::cout << std::endl;