pktools 2.6.7
Processing Kernel for geospatial data
pkfssvm.cc
1/**********************************************************************
2pkfssvm.cc: feature selection for support vector machine classifier pksvm
3Copyright (C) 2008-2014 Pieter Kempeneers
4
5This file is part of pktools
6
7pktools is free software: you can redistribute it and/or modify
8it under the terms of the GNU General Public License as published by
9the Free Software Foundation, either version 3 of the License, or
10(at your option) any later version.
11
12pktools is distributed in the hope that it will be useful,
13but WITHOUT ANY WARRANTY; without even the implied warranty of
14MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15GNU General Public License for more details.
16
17You should have received a copy of the GNU General Public License
18along with pktools. If not, see <http://www.gnu.org/licenses/>.
19***********************************************************************/
20#include <stdlib.h>
21#include <vector>
22#include <string>
23#include <map>
24#include <algorithm>
25#include "base/Optionpk.h"
26#include "algorithms/ConfusionMatrix.h"
27#include "algorithms/CostFactorySVM.h"
28#include "algorithms/FeatureSelector.h"
29#include "algorithms/svm.h"
30#include "imageclasses/ImgReaderOgr.h"
31
32#ifdef HAVE_CONFIG_H
33#include <config.h>
34#endif
35
36/******************************************************************************/
98using namespace std;
99
100enum SelectorValue { NA=0, SFFS=1, SFS=2, SBS=3, BFS=4};
101
102// CostFactorySVM::CostFactorySVM()
103// : CostFactory(2,0), m_svm_type("C_SVC"), m_kernel_type("radial"), m_kernel_degree(3), m_gamma(1.0), m_coef0(0), m_ccost(1000), m_nu(0.5), m_epsilon_loss(100), m_cache(100), m_epsilon_tol(0.001), m_shrinking(false), m_prob_est(true){
104// }
105
106// CostFactorySVM::~CostFactorySVM(){
107// }
108
109// CostFactorySVM::CostFactorySVM(std::string svm_type, std::string kernel_type, unsigned short kernel_degree, float gamma, float coef0, float ccost, float nu, float epsilon_loss, int cache, float epsilon_tol, bool shrinking, bool prob_est, unsigned short cv, bool verbose)
110// : CostFactory(cv,verbose), m_svm_type(svm_type), m_kernel_type(kernel_type), m_kernel_degree(kernel_degree), m_gamma(gamma), m_coef0(coef0), m_ccost(ccost), m_nu(nu), m_epsilon_loss(epsilon_loss), m_cache(cache), m_epsilon_tol(epsilon_tol), m_shrinking(shrinking), m_prob_est(prob_est){};
111
112// double CostFactorySVM::getCost(const vector<Vector2d<float> > &trainingFeatures){
113// std::map<std::string, svm::SVM_TYPE> svmMap;
114
115// svmMap["C_SVC"]=svm::C_SVC;
116// svmMap["nu_SVC"]=svm::nu_SVC;
117// svmMap["one_class"]=svm::one_class;
118// svmMap["epsilon_SVR"]=svm::epsilon_SVR;
119// svmMap["nu_SVR"]=svm::nu_SVR;
120
121// std::map<std::string, svm::KERNEL_TYPE> kernelMap;
122
123// kernelMap["linear"]=svm::linear;
124// kernelMap["polynomial"]=svm::polynomial;
125// kernelMap["radial"]=svm::radial;
126// kernelMap["sigmoid;"]=svm::sigmoid;
127
128// unsigned short nclass=trainingFeatures.size();
129// unsigned int ntraining=0;
130// unsigned int ntest=0;
131// for(int iclass=0;iclass<nclass;++iclass){
132// ntraining+=m_nctraining[iclass];
133// ntest+=m_nctest[iclass];
134// }
135// if(ntest)
136// assert(!m_cv);
137// if(!m_cv)
138// assert(ntest);
139// unsigned short nFeatures=trainingFeatures[0][0].size();
140
141// struct svm_parameter param;
142// param.svm_type = svmMap[m_svm_type];
143// param.kernel_type = kernelMap[m_kernel_type];
144// param.degree = m_kernel_degree;
145// param.gamma = (m_gamma>0)? m_gamma : 1.0/nFeatures;
146// param.coef0 = m_coef0;
147// param.nu = m_nu;
148// param.cache_size = m_cache;
149// param.C = m_ccost;
150// param.eps = m_epsilon_tol;
151// param.p = m_epsilon_loss;
152// param.shrinking = (m_shrinking)? 1 : 0;
153// param.probability = (m_prob_est)? 1 : 0;
154// param.nr_weight = 0;//not used: I use priors and balancing
155// param.weight_label = NULL;
156// param.weight = NULL;
157// param.verbose=(m_verbose>1)? true:false;
158// struct svm_model* svm;
159// struct svm_problem prob;
160// struct svm_node* x_space;
161
162// prob.l=ntraining;
163// prob.y = Malloc(double,prob.l);
164// prob.x = Malloc(struct svm_node *,prob.l);
165// x_space = Malloc(struct svm_node,(nFeatures+1)*ntraining);
166// unsigned long int spaceIndex=0;
167// int lIndex=0;
168// for(int iclass=0;iclass<nclass;++iclass){
169// // for(int isample=0;isample<trainingFeatures[iclass].size();++isample){
170// for(int isample=0;isample<m_nctraining[iclass];++isample){
171// prob.x[lIndex]=&(x_space[spaceIndex]);
172// for(int ifeature=0;ifeature<nFeatures;++ifeature){
173// x_space[spaceIndex].index=ifeature+1;
174// x_space[spaceIndex].value=trainingFeatures[iclass][isample][ifeature];
175// ++spaceIndex;
176// }
177// x_space[spaceIndex++].index=-1;
178// prob.y[lIndex]=iclass;
179// ++lIndex;
180// }
181// }
182
183// assert(lIndex==prob.l);
184// if(m_verbose>2)
185// std::cout << "checking parameters" << std::endl;
186// svm_check_parameter(&prob,&param);
187// if(m_verbose>2)
188// std::cout << "parameters ok, training" << std::endl;
189// svm=svm_train(&prob,&param);
190// if(m_verbose>2)
191// std::cout << "SVM is now trained" << std::endl;
192
193// m_cm.clearResults();
194// if(m_cv>1){
195// double *target = Malloc(double,prob.l);
196// svm_cross_validation(&prob,&param,m_cv,target);
197// assert(param.svm_type != EPSILON_SVR&&param.svm_type != NU_SVR);//only for regression
198// for(int i=0;i<prob.l;i++){
199// string refClassName=m_nameVector[prob.y[i]];
200// string className=m_nameVector[target[i]];
201// if(m_classValueMap.size())
202// m_cm.incrementResult(type2string<short>(m_classValueMap[refClassName]),type2string<short>(m_classValueMap[className]),1.0);
203// else
204// m_cm.incrementResult(m_cm.getClass(prob.y[i]),m_cm.getClass(target[i]),1.0);
205// }
206// free(target);
207// }
208// else{
209// struct svm_node *x_test;
210// vector<double> result(nclass);
211// x_test = Malloc(struct svm_node,(nFeatures+1));
212// for(int iclass=0;iclass<nclass;++iclass){
213// for(int isample=0;isample<m_nctest[iclass];++isample){
214// for(int ifeature=0;ifeature<nFeatures;++ifeature){
215// x_test[ifeature].index=ifeature+1;
216// x_test[ifeature].value=trainingFeatures[iclass][m_nctraining[iclass]+isample][ifeature];
217// }
218// x_test[nFeatures].index=-1;
219// double predict_label=0;
220// assert(svm_check_probability_model(svm));
221// predict_label = svm_predict_probability(svm,x_test,&(result[0]));
222// // predict_label = svm_predict(svm,x_test);
223// string refClassName=m_nameVector[iclass];
224// string className=m_nameVector[static_cast<short>(predict_label)];
225// if(m_classValueMap.size())
226// m_cm.incrementResult(type2string<short>(m_classValueMap[refClassName]),type2string<short>(m_classValueMap[className]),1.0);
227// else
228// m_cm.incrementResult(refClassName,className,1.0);
229// }
230// }
231// free(x_test);
232// }
233// if(m_verbose>1)
234// std::cout << m_cm << std::endl;
235// assert(m_cm.nReference());
236// // if(m_verbose)
237
238// // std::cout << m_cm << std::endl;
239// // std::cout << "Kappa: " << m_cm.kappa() << std::endl;
240// // double se95_oa=0;
241// // double doa=0;
242// // doa=m_cm.oa_pct(&se95_oa);
243// // std::cout << "Overall Accuracy: " << doa << " (" << se95_oa << ")" << std::endl;
244
245// // *NOTE* Because svm_model contains pointers to svm_problem, you can
246// // not free the memory used by svm_problem if you are still using the
247// // svm_model produced by svm_train().
248// // however, we will re-train the svm later on after the feature selection
249// free(prob.y);
250// free(prob.x);
251// free(x_space);
252// svm_free_and_destroy_model(&(svm));
253
254// return(m_cm.kappa());
255// }
256
257int main(int argc, char *argv[])
258{
259 // vector<double> priors;
260
261 //--------------------------- command line options ------------------------------------
262 Optionpk<string> input_opt("i", "input", "input test set (leave empty to perform a cross validation based on training only)");
263 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).");
264 Optionpk<string> tlayer_opt("tln", "tln", "training layer name(s)");
265 Optionpk<string> label_opt("label", "label", "identifier for class label in training vector file.","label");
266 Optionpk<unsigned short> maxFeatures_opt("n", "nf", "number of features to select (0 to select optimal number, see also ecost option)", 0);
267 Optionpk<unsigned int> balance_opt("bal", "balance", "balance the input data to this number of samples for each class", 0);
268 Optionpk<bool> random_opt("random","random", "in case of balance, randomize input data", true);
269 Optionpk<int> minSize_opt("min", "min", "if number of training pixels is less then min, do not take this class into account", 0);
270 Optionpk<unsigned short> band_opt("b", "band", "band index (starting from 0, either use band option or use start to end)");
271 Optionpk<unsigned short> bstart_opt("sband", "startband", "Start band sequence number");
272 Optionpk<unsigned short> bend_opt("eband", "endband", "End band sequence number");
273 Optionpk<double> offset_opt("offset", "offset", "offset value for each spectral band input features: refl[band]=(DN[band]-offset[band])/scale[band]", 0.0);
274 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);
275 Optionpk<string> selector_opt("sm", "sm", "feature selection method (sffs=sequential floating forward search,sfs=sequential forward search, sbs, sequential backward search ,bfs=brute force search)","sffs");
276 Optionpk<float> epsilon_cost_opt("ecost", "ecost", "epsilon for stopping criterion in cost function to determine optimal number of features",0.001);
277
278 Optionpk<std::string> svm_type_opt("svmt", "svmtype", "type of SVM (C_SVC, nu_SVC,one_class, epsilon_SVR, nu_SVR)","C_SVC");
279 Optionpk<std::string> kernel_type_opt("kt", "kerneltype", "type of kernel function (linear,polynomial,radial,sigmoid) ","radial");
280 Optionpk<unsigned short> kernel_degree_opt("kd", "kd", "degree in kernel function",3);
281 Optionpk<float> gamma_opt("g", "gamma", "gamma in kernel function",1.0);
282 Optionpk<float> coef0_opt("c0", "coef0", "coef0 in kernel function",0);
283 Optionpk<float> ccost_opt("cc", "ccost", "the parameter C of C-SVC, epsilon-SVR, and nu-SVR",1000);
284 Optionpk<float> nu_opt("nu", "nu", "the parameter nu of nu-SVC, one-class SVM, and nu-SVR",0.5);
285 Optionpk<float> epsilon_loss_opt("eloss", "eloss", "the epsilon in loss function of epsilon-SVR",0.1);
286 Optionpk<int> cache_opt("cache", "cache", "cache memory size in MB",100);
287 Optionpk<float> epsilon_tol_opt("etol", "etol", "the tolerance of termination criterion",0.001);
288 Optionpk<bool> shrinking_opt("shrink", "shrink", "whether to use the shrinking heuristics",false);
289 Optionpk<bool> prob_est_opt("pe", "probest", "whether to train a SVC or SVR model for probability estimates",true,2);
290 Optionpk<unsigned short> cv_opt("cv", "cv", "n-fold cross validation mode",2);
291 Optionpk<string> classname_opt("c", "class", "list of class names.");
292 Optionpk<short> classvalue_opt("r", "reclass", "list of class values (use same order as in classname opt.");
293 Optionpk<short> verbose_opt("v", "verbose", "set to: 0 (results only), 1 (confusion matrix), 2 (debug)",0,2);
294
295 tlayer_opt.setHide(1);
296 label_opt.setHide(1);
297 balance_opt.setHide(1);
298 random_opt.setHide(1);
299 minSize_opt.setHide(1);
300 band_opt.setHide(1);
301 bstart_opt.setHide(1);
302 bend_opt.setHide(1);
303 offset_opt.setHide(1);
304 scale_opt.setHide(1);
305 svm_type_opt.setHide(1);
306 kernel_type_opt.setHide(1);
307 kernel_degree_opt.setHide(1);
308 gamma_opt.setHide(1);
309 coef0_opt.setHide(1);
310 ccost_opt.setHide(1);
311 nu_opt.setHide(1);
312 epsilon_loss_opt.setHide(1);
313 cache_opt.setHide(1);
314 epsilon_tol_opt.setHide(1);
315 shrinking_opt.setHide(1);
316 prob_est_opt.setHide(1);
317 selector_opt.setHide(1);
318 epsilon_cost_opt.setHide(1);
319 cv_opt.setHide(1);
320 classname_opt.setHide(1);
321 classvalue_opt.setHide(1);
322
323 bool doProcess;//stop process when program was invoked with help option (-h --help)
324 try{
325 doProcess=input_opt.retrieveOption(argc,argv);
326 training_opt.retrieveOption(argc,argv);
327 maxFeatures_opt.retrieveOption(argc,argv);
328 tlayer_opt.retrieveOption(argc,argv);
329 label_opt.retrieveOption(argc,argv);
330 balance_opt.retrieveOption(argc,argv);
331 random_opt.retrieveOption(argc,argv);
332 minSize_opt.retrieveOption(argc,argv);
333 band_opt.retrieveOption(argc,argv);
334 bstart_opt.retrieveOption(argc,argv);
335 bend_opt.retrieveOption(argc,argv);
336 offset_opt.retrieveOption(argc,argv);
337 scale_opt.retrieveOption(argc,argv);
338 svm_type_opt.retrieveOption(argc,argv);
339 kernel_type_opt.retrieveOption(argc,argv);
340 kernel_degree_opt.retrieveOption(argc,argv);
341 gamma_opt.retrieveOption(argc,argv);
342 coef0_opt.retrieveOption(argc,argv);
343 ccost_opt.retrieveOption(argc,argv);
344 nu_opt.retrieveOption(argc,argv);
345 epsilon_loss_opt.retrieveOption(argc,argv);
346 cache_opt.retrieveOption(argc,argv);
347 epsilon_tol_opt.retrieveOption(argc,argv);
348 shrinking_opt.retrieveOption(argc,argv);
349 prob_est_opt.retrieveOption(argc,argv);
350 selector_opt.retrieveOption(argc,argv);
351 epsilon_cost_opt.retrieveOption(argc,argv);
352 cv_opt.retrieveOption(argc,argv);
353 classname_opt.retrieveOption(argc,argv);
354 classvalue_opt.retrieveOption(argc,argv);
355 verbose_opt.retrieveOption(argc,argv);
356 }
357 catch(string predefinedString){
358 std::cout << predefinedString << std::endl;
359 exit(0);
360 }
361 if(!doProcess){
362 cout << endl;
363 cout << "Usage: pkfssvm -t training -n number" << endl;
364 cout << endl;
365 std::cout << "short option -h shows basic options only, use long option --help to show all options" << std::endl;
366 exit(0);//help was invoked, stop processing
367 }
368
369 CostFactorySVM costfactory(svm_type_opt[0], kernel_type_opt[0], kernel_degree_opt[0], gamma_opt[0], coef0_opt[0], ccost_opt[0], 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]);
370
371 assert(training_opt.size());
372 if(input_opt.size())
373 costfactory.setCv(0);
374 if(verbose_opt[0]>=1){
375 if(input_opt.size())
376 std::cout << "input filename: " << input_opt[0] << std::endl;
377 std::cout << "training vector file: " << std::endl;
378 for(int ifile=0;ifile<training_opt.size();++ifile)
379 std::cout << training_opt[ifile] << std::endl;
380 std::cout << "verbose: " << verbose_opt[0] << std::endl;
381 }
382
383 static std::map<std::string, SelectorValue> selMap;
384 //initialize selMap
385 selMap["sffs"]=SFFS;
386 selMap["sfs"]=SFS;
387 selMap["sbs"]=SBS;
388 selMap["bfs"]=BFS;
389
390 unsigned int totalSamples=0;
391 unsigned int totalTestSamples=0;
392
393 unsigned short nclass=0;
394 int nband=0;
395 int startBand=2;//first two bands represent X and Y pos
396
397 // if(priors_opt.size()>1){//priors from argument list
398 // priors.resize(priors_opt.size());
399 // double normPrior=0;
400 // for(int iclass=0;iclass<priors_opt.size();++iclass){
401 // priors[iclass]=priors_opt[iclass];
402 // normPrior+=priors[iclass];
403 // }
404 // //normalize
405 // for(int iclass=0;iclass<priors_opt.size();++iclass)
406 // priors[iclass]/=normPrior;
407 // }
408
409 //convert start and end band options to vector of band indexes
410 try{
411 if(bstart_opt.size()){
412 if(bend_opt.size()!=bstart_opt.size()){
413 string errorstring="Error: options for start and end band indexes must be provided as pairs, missing end band";
414 throw(errorstring);
415 }
416 band_opt.clear();
417 for(int ipair=0;ipair<bstart_opt.size();++ipair){
418 if(bend_opt[ipair]<=bstart_opt[ipair]){
419 string errorstring="Error: index for end band must be smaller then start band";
420 throw(errorstring);
421 }
422 for(int iband=bstart_opt[ipair];iband<=bend_opt[ipair];++iband)
423 band_opt.push_back(iband);
424 }
425 }
426 }
427 catch(string error){
428 cerr << error << std::endl;
429 exit(1);
430 }
431 //sort bands
432 if(band_opt.size())
433 std::sort(band_opt.begin(),band_opt.end());
434
435 if(classname_opt.size()){
436 assert(classname_opt.size()==classvalue_opt.size());
437 for(int iclass=0;iclass<classname_opt.size();++iclass)
438 costfactory.setClassValueMap(classname_opt[iclass],classvalue_opt[iclass]);
439 }
440
441 //----------------------------------- Training -------------------------------
442 vector<double> offset;
443 vector<double> scale;
444 vector< Vector2d<float> > trainingPixels;//[class][sample][band]
445 vector< Vector2d<float> > testPixels;//[class][sample][band]
446 map<string,Vector2d<float> > trainingMap;
447 map<string,Vector2d<float> > testMap;
448 vector<string> fields;
449
450 struct svm_problem prob;
451 //organize training data
452 trainingPixels.clear();
453 testPixels.clear();
454 if(verbose_opt[0]>=1)
455 std::cout << "reading training file " << training_opt[0] << std::endl;
456 try{
457 ImgReaderOgr trainingReader(training_opt[0]);
458 if(band_opt.size()){
459 totalSamples=trainingReader.readDataImageOgr(trainingMap,fields,band_opt,label_opt[0],tlayer_opt,verbose_opt[0]);
460 if(input_opt.size()){
461 ImgReaderOgr inputReader(input_opt[0]);
462 totalTestSamples=inputReader.readDataImageOgr(testMap,fields,band_opt,label_opt[0],tlayer_opt,verbose_opt[0]);
463 inputReader.close();
464 }
465 }
466 else{
467 totalSamples=trainingReader.readDataImageOgr(trainingMap,fields,0,0,label_opt[0],tlayer_opt,verbose_opt[0]);
468 if(input_opt.size()){
469 ImgReaderOgr inputReader(input_opt[0]);
470 totalTestSamples=inputReader.readDataImageOgr(testMap,fields,0,0,label_opt[0],tlayer_opt,verbose_opt[0]);
471 inputReader.close();
472 }
473 }
474 if(trainingMap.size()<2){
475 string errorstring="Error: could not read at least two classes from training input file";
476 throw(errorstring);
477 }
478 if(input_opt.size()&&testMap.size()<2){
479 string errorstring="Error: could not read at least two classes from test input file";
480 throw(errorstring);
481 }
482 trainingReader.close();
483 }
484 catch(string error){
485 cerr << error << std::endl;
486 exit(1);
487 }
488 catch(std::exception& e){
489 std::cerr << "Error: ";
490 std::cerr << e.what() << std::endl;
491 std::cerr << CPLGetLastErrorMsg() << std::endl;
492 exit(1);
493 }
494 catch(...){
495 cerr << "error caught" << std::endl;
496 exit(1);
497 }
498 //todo: delete class 0 ?
499 // if(verbose_opt[0]>=1)
500 // std::cout << "erasing class 0 from training set (" << trainingMap[0].size() << " from " << totalSamples << ") samples" << std::endl;
501 // totalSamples-=trainingMap[0].size();
502 // trainingMap.erase(0);
503
504 if(verbose_opt[0]>1)
505 std::cout << "training pixels: " << std::endl;
506 map<string,Vector2d<float> >::iterator mapit=trainingMap.begin();
507 while(mapit!=trainingMap.end()){
508 //delete small classes
509 if((mapit->second).size()<minSize_opt[0]){
510 trainingMap.erase(mapit);
511 continue;
512 }
513 costfactory.pushBackName(mapit->first);
514 trainingPixels.push_back(mapit->second);
515 if(verbose_opt[0]>1)
516 std::cout << mapit->first << ": " << (mapit->second).size() << " samples" << std::endl;
517 ++mapit;
518 }
519 nclass=trainingPixels.size();
520 if(classname_opt.size())
521 assert(nclass==classname_opt.size());
522 nband=trainingPixels[0][0].size()-2;//X and Y//trainingPixels[0][0].size();
523
524 mapit=testMap.begin();
525 while(mapit!=testMap.end()){
526 if(costfactory.getClassValueMap().size()){
527 // if(classValueMap.size()){
528 //check if name in test is covered by classname_opt (values can not be 0)
529 if((costfactory.getClassValueMap())[mapit->first]>0){
530 ;//ok, no need to print to std::cout
531 }
532 else{
533 std::cerr << "Error: names in classname option are not complete, please check names in test vector and make sure classvalue is > 0" << std::endl;
534 exit(1);
535 }
536 }
537 //no need to delete small classes for test sample
538 testPixels.push_back(mapit->second);
539 if(verbose_opt[0]>1)
540 std::cout << mapit->first << ": " << (mapit->second).size() << " samples" << std::endl;
541 ++mapit;
542 }
543 if(input_opt.size()){
544 assert(nclass==testPixels.size());
545 assert(nband=testPixels[0][0].size()-2);//X and Y//testPixels[0][0].size();
546 assert(!cv_opt[0]);
547 }
548
549 //do not remove outliers here: could easily be obtained through ogr2ogr -where 'B2<110' output.shp input.shp
550 //balance training data
551 //todo: do I mean to use random_opt?
552 if(balance_opt[0]>0){
553 if(random_opt[0])
554 srand(time(NULL));
555 totalSamples=0;
556 for(int iclass=0;iclass<nclass;++iclass){
557 if(trainingPixels[iclass].size()>balance_opt[0]){
558 while(trainingPixels[iclass].size()>balance_opt[0]){
559 int index=rand()%trainingPixels[iclass].size();
560 trainingPixels[iclass].erase(trainingPixels[iclass].begin()+index);
561 }
562 }
563 else{
564 int oldsize=trainingPixels[iclass].size();
565 for(int isample=trainingPixels[iclass].size();isample<balance_opt[0];++isample){
566 int index = rand()%oldsize;
567 trainingPixels[iclass].push_back(trainingPixels[iclass][index]);
568 }
569 }
570 totalSamples+=trainingPixels[iclass].size();
571 }
572 assert(totalSamples==nclass*balance_opt[0]);
573 }
574
575 //set scale and offset
576 offset.resize(nband);
577 scale.resize(nband);
578 if(offset_opt.size()>1)
579 assert(offset_opt.size()==nband);
580 if(scale_opt.size()>1)
581 assert(scale_opt.size()==nband);
582 for(int iband=0;iband<nband;++iband){
583 if(verbose_opt[0]>1)
584 std::cout << "scaling for band" << iband << std::endl;
585 offset[iband]=(offset_opt.size()==1)?offset_opt[0]:offset_opt[iband];
586 scale[iband]=(scale_opt.size()==1)?scale_opt[0]:scale_opt[iband];
587 //search for min and maximum
588 if(scale[iband]<=0){
589 float theMin=trainingPixels[0][0][iband+startBand];
590 float theMax=trainingPixels[0][0][iband+startBand];
591 for(int iclass=0;iclass<nclass;++iclass){
592 for(int isample=0;isample<trainingPixels[iclass].size();++isample){
593 if(theMin>trainingPixels[iclass][isample][iband+startBand])
594 theMin=trainingPixels[iclass][isample][iband+startBand];
595 if(theMax<trainingPixels[iclass][isample][iband+startBand])
596 theMax=trainingPixels[iclass][isample][iband+startBand];
597 }
598 }
599 offset[iband]=theMin+(theMax-theMin)/2.0;
600 scale[iband]=(theMax-theMin)/2.0;
601 if(verbose_opt[0]>1){
602 std::cout << "Extreme image values for band " << iband << ": [" << theMin << "," << theMax << "]" << std::endl;
603 std::cout << "Using offset, scale: " << offset[iband] << ", " << scale[iband] << std::endl;
604 std::cout << "scaled values for band " << iband << ": [" << (theMin-offset[iband])/scale[iband] << "," << (theMax-offset[iband])/scale[iband] << "]" << std::endl;
605 }
606 }
607 }
608
609 // if(priors_opt.size()==1){//default: equal priors for each class
610 // priors.resize(nclass);
611 // for(int iclass=0;iclass<nclass;++iclass)
612 // priors[iclass]=1.0/nclass;
613 // }
614 // assert(priors_opt.size()==1||priors_opt.size()==nclass);
615
616 if(verbose_opt[0]>=1){
617 std::cout << "number of bands: " << nband << std::endl;
618 std::cout << "number of classes: " << nclass << std::endl;
619 // std::cout << "priors:";
620 // for(int iclass=0;iclass<nclass;++iclass)
621 // std::cout << " " << priors[iclass];
622 // std::cout << std::endl;
623 }
624
625 //set names in confusion matrix using nameVector
626 vector<string> nameVector=costfactory.getNameVector();
627 for(int iname=0;iname<nameVector.size();++iname){
628 if(costfactory.getClassValueMap().empty())
629 costfactory.pushBackClassName(nameVector[iname]);
630 // cm.pushBackClassName(nameVector[iname]);
631 else if(costfactory.getClassIndex(type2string<short>((costfactory.getClassValueMap())[nameVector[iname]]))<0)
632 costfactory.pushBackClassName(type2string<short>((costfactory.getClassValueMap())[nameVector[iname]]));
633 }
634
635
636 //Calculate features of training (and test) set
637
638 vector<unsigned int> nctraining;
639 vector<unsigned int> nctest;
640 nctraining.resize(nclass);
641 nctest.resize(nclass);
642 vector< Vector2d<float> > trainingFeatures(nclass);
643 for(int iclass=0;iclass<nclass;++iclass){
644 if(verbose_opt[0]>=1)
645 std::cout << "calculating features for class " << iclass << std::endl;
646 nctraining[iclass]=trainingPixels[iclass].size();
647 if(verbose_opt[0]>=1)
648 std::cout << "nctraining[" << iclass << "]: " << nctraining[iclass] << std::endl;
649 if(testPixels.size()>iclass){
650 nctest[iclass]=testPixels[iclass].size();
651 if(verbose_opt[0]>=1){
652 std::cout << "nctest[" << iclass << "]: " << nctest[iclass] << std::endl;
653 }
654 }
655 else
656 nctest[iclass]=0;
657
658 trainingFeatures[iclass].resize(nctraining[iclass]+nctest[iclass]);
659 for(int isample=0;isample<nctraining[iclass];++isample){
660 //scale pixel values according to scale and offset!!!
661 for(int iband=0;iband<nband;++iband){
662 assert(trainingPixels[iclass].size()>isample);
663 assert(trainingPixels[iclass][isample].size()>iband+startBand);
664 assert(offset.size()>iband);
665 assert(scale.size()>iband);
666 float value=trainingPixels[iclass][isample][iband+startBand];
667 trainingFeatures[iclass][isample].push_back((value-offset[iband])/scale[iband]);
668 }
669 }
670 for(int isample=0;isample<nctest[iclass];++isample){
671 //scale pixel values according to scale and offset!!!
672 for(int iband=0;iband<nband;++iband){
673 assert(testPixels[iclass].size()>isample);
674 assert(testPixels[iclass][isample].size()>iband+startBand);
675 assert(offset.size()>iband);
676 assert(scale.size()>iband);
677 float value=testPixels[iclass][isample][iband+startBand];
678 // testFeatures[iclass][isample].push_back((value-offset[iband])/scale[iband]);
679 trainingFeatures[iclass][nctraining[iclass]+isample].push_back((value-offset[iband])/scale[iband]);
680 }
681 }
682 assert(trainingFeatures[iclass].size()==nctraining[iclass]+nctest[iclass]);
683 }
684
685 costfactory.setNcTraining(nctraining);
686 costfactory.setNcTest(nctest);
687 int nFeatures=trainingFeatures[0][0].size();
688 int maxFeatures=(maxFeatures_opt[0])? maxFeatures_opt[0] : 1;
689 double previousCost=-1;
690 double cost=0;
691 list<int> subset;//set of selected features (levels) for each class combination
692 FeatureSelector selector;
693 try{
694 if(maxFeatures>=nFeatures){
695 subset.clear();
696 for(int ifeature=0;ifeature<nFeatures;++ifeature)
697 subset.push_back(ifeature);
698 cost=costfactory.getCost(trainingFeatures);
699 }
700 else{
701 while(fabs(cost-previousCost)>=epsilon_cost_opt[0]){
702 previousCost=cost;
703 switch(selMap[selector_opt[0]]){
704 case(SFFS):
705 subset.clear();//needed to clear in case of floating and brute force search
706 cost=selector.floating(trainingFeatures,costfactory,subset,maxFeatures,epsilon_cost_opt[0],verbose_opt[0]);
707 break;
708 case(SFS):
709 cost=selector.forward(trainingFeatures,costfactory,subset,maxFeatures,verbose_opt[0]);
710 break;
711 case(SBS):
712 cost=selector.backward(trainingFeatures,costfactory,subset,maxFeatures,verbose_opt[0]);
713 break;
714 case(BFS):
715 subset.clear();//needed to clear in case of floating and brute force search
716 cost=selector.bruteForce(trainingFeatures,costfactory,subset,maxFeatures,verbose_opt[0]);
717 break;
718 default:
719 std::cout << "Error: selector not supported, please use sffs, sfs, sbs or bfs" << std::endl;
720 exit(1);
721 break;
722 }
723 if(verbose_opt[0]>1){
724 std::cout << "cost: " << cost << std::endl;
725 std::cout << "previousCost: " << previousCost << std::endl;
726 std::cout << std::setprecision(12) << "cost-previousCost: " << cost - previousCost << " ( " << epsilon_cost_opt[0] << ")" << std::endl;
727 }
728 if(!maxFeatures_opt[0])
729 ++maxFeatures;
730 else
731 break;
732 }
733 }
734 }
735 catch(...){
736 std::cout << "caught feature selection" << std::endl;
737 exit(1);
738 }
739
740 if(verbose_opt[0])
741 cout <<"cost: " << cost << endl;
742 subset.sort();
743 for(list<int>::const_iterator lit=subset.begin();lit!=subset.end();++lit)
744 std::cout << " -b " << *lit;
745 std::cout << std::endl;
746 // if((*(lit))!=subset.back())
747 // else
748 // cout << endl;
749
750 // *NOTE* Because svm_model contains pointers to svm_problem, you can
751 // not free the memory used by svm_problem if you are still using the
752 // svm_model produced by svm_train().
753
754 // free(prob.y);
755 // free(prob.x);
756 // free(x_space);
757 // svm_destroy_param(&param);
758 return 0;
759}
760