24 #include "imageclasses/ImgReaderGdal.h"
25 #include "imageclasses/ImgWriterGdal.h"
26 #include "imageclasses/ImgReaderOgr.h"
27 #include "imageclasses/ImgWriterOgr.h"
28 #include "base/Optionpk.h"
29 #include "base/PosValue.h"
30 #include "algorithms/ConfusionMatrix.h"
31 #include "algorithms/svm.h"
119 enum SVM_TYPE {C_SVC=0, nu_SVC=1,one_class=2, epsilon_SVR=3, nu_SVR=4};
120 enum KERNEL_TYPE {linear=0,polynomial=1,radial=2,sigmoid=3};
123 #define Malloc(type,n) (type *)malloc((n)*sizeof(type))
127 int main(
int argc,
char *argv[])
129 vector<double> priors;
133 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). Use multiple training files for bootstrap aggregation (alternative to the bag and bsize options, where a random subset is taken from a single training file)");
135 Optionpk<string> label_opt(
"label",
"label",
"Attribute name for class label in training vector file.",
"label");
136 Optionpk<unsigned int> balance_opt(
"bal",
"balance",
"Balance the input data to this number of samples for each class", 0);
137 Optionpk<bool> random_opt(
"random",
"random",
"Randomize training data for balancing and bagging",
true, 2);
138 Optionpk<int> minSize_opt(
"min",
"min",
"If number of training pixels is less then min, do not take this class into account (0: consider all classes)", 0);
139 Optionpk<unsigned short> band_opt(
"b",
"band",
"Band index (starting from 0, either use band option or use start to end)");
142 Optionpk<double> offset_opt(
"offset",
"offset",
"Offset value for each spectral band input features: refl[band]=(DN[band]-offset[band])/scale[band]", 0.0);
143 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);
144 Optionpk<double> priors_opt(
"prior",
"prior",
"Prior probabilities for each class (e.g., -p 0.3 -p 0.3 -p 0.2 ). Used for input only (ignored for cross validation)", 0.0);
145 Optionpk<string> priorimg_opt(
"pim",
"priorimg",
"Prior probability image (multi-band img with band for each class",
"",2);
147 Optionpk<string> cmformat_opt(
"cmf",
"cmf",
"Format for confusion matrix (ascii or latex)",
"ascii");
148 Optionpk<std::string> svm_type_opt(
"svmt",
"svmtype",
"Type of SVM (C_SVC, nu_SVC,one_class, epsilon_SVR, nu_SVR)",
"C_SVC");
149 Optionpk<std::string> kernel_type_opt(
"kt",
"kerneltype",
"Type of kernel function (linear,polynomial,radial,sigmoid) ",
"radial");
151 Optionpk<float> gamma_opt(
"g",
"gamma",
"Gamma in kernel function",1.0);
152 Optionpk<float> coef0_opt(
"c0",
"coef0",
"Coef0 in kernel function",0);
153 Optionpk<float> ccost_opt(
"cc",
"ccost",
"The parameter C of C_SVC, epsilon_SVR, and nu_SVR",1000);
154 Optionpk<float> nu_opt(
"nu",
"nu",
"The parameter nu of nu_SVC, one_class SVM, and nu_SVR",0.5);
155 Optionpk<float> epsilon_loss_opt(
"eloss",
"eloss",
"The epsilon in loss function of epsilon_SVR",0.1);
156 Optionpk<int> cache_opt(
"cache",
"cache",
"Cache memory size in MB",100);
157 Optionpk<float> epsilon_tol_opt(
"etol",
"etol",
"The tolerance of termination criterion",0.001);
158 Optionpk<bool> shrinking_opt(
"shrink",
"shrink",
"Whether to use the shrinking heuristics",
false);
159 Optionpk<bool> prob_est_opt(
"pe",
"probest",
"Whether to train a SVC or SVR model for probability estimates",
true,2);
161 Optionpk<unsigned short> comb_opt(
"comb",
"comb",
"How to combine bootstrap aggregation classifiers (0: sum rule, 1: product rule, 2: max rule). Also used to aggregate classes with rc option.",0);
163 Optionpk<int> bagSize_opt(
"bagsize",
"bagsize",
"Percentage of features used from available training features for each bootstrap aggregation (one size for all classes, or a different size for each class respectively", 100);
164 Optionpk<string> classBag_opt(
"cb",
"classbag",
"Output for each individual bootstrap aggregation");
165 Optionpk<string> mask_opt(
"m",
"mask",
"Only classify within specified mask (vector or raster). For raster mask, set nodata values with the option msknodata.");
166 Optionpk<short> msknodata_opt(
"msknodata",
"msknodata",
"Mask value(s) not to consider for classification. Values will be taken over in classification image.", 0);
169 Optionpk<string> oformat_opt(
"of",
"oformat",
"Output image format (see also gdal_translate).",
"GTiff");
170 Optionpk<string> option_opt(
"co",
"co",
"Creation option for output file. Multiple options can be specified.");
171 Optionpk<string> colorTable_opt(
"ct",
"ct",
"Color table in ASCII format having 5 columns: id R G B ALFA (0: transparent, 255: solid)");
173 Optionpk<string> entropy_opt(
"entropy",
"entropy",
"Entropy image (measure for uncertainty of classifier output",
"",2);
174 Optionpk<string> active_opt(
"active",
"active",
"Ogr output for active training sample.",
"",2);
175 Optionpk<string> ogrformat_opt(
"f",
"f",
"Output ogr format for active training sample",
"SQLite");
178 Optionpk<short> classvalue_opt(
"r",
"reclass",
"List of class values (use same order as in class opt).");
181 oformat_opt.setHide(1);
182 option_opt.setHide(1);
184 bstart_opt.setHide(1);
186 balance_opt.setHide(1);
187 minSize_opt.setHide(1);
189 bagSize_opt.setHide(1);
191 classBag_opt.setHide(1);
193 priorimg_opt.setHide(1);
194 offset_opt.setHide(1);
195 scale_opt.setHide(1);
196 svm_type_opt.setHide(1);
197 kernel_type_opt.setHide(1);
198 kernel_degree_opt.setHide(1);
199 coef0_opt.setHide(1);
201 epsilon_loss_opt.setHide(1);
202 cache_opt.setHide(1);
203 epsilon_tol_opt.setHide(1);
204 shrinking_opt.setHide(1);
205 prob_est_opt.setHide(1);
206 entropy_opt.setHide(1);
207 active_opt.setHide(1);
208 nactive_opt.setHide(1);
209 random_opt.setHide(1);
211 verbose_opt.setHide(2);
215 doProcess=training_opt.retrieveOption(argc,argv);
216 input_opt.retrieveOption(argc,argv);
217 output_opt.retrieveOption(argc,argv);
218 cv_opt.retrieveOption(argc,argv);
219 cmformat_opt.retrieveOption(argc,argv);
220 tlayer_opt.retrieveOption(argc,argv);
221 classname_opt.retrieveOption(argc,argv);
222 classvalue_opt.retrieveOption(argc,argv);
223 oformat_opt.retrieveOption(argc,argv);
224 ogrformat_opt.retrieveOption(argc,argv);
225 option_opt.retrieveOption(argc,argv);
226 colorTable_opt.retrieveOption(argc,argv);
227 label_opt.retrieveOption(argc,argv);
228 priors_opt.retrieveOption(argc,argv);
229 gamma_opt.retrieveOption(argc,argv);
230 ccost_opt.retrieveOption(argc,argv);
231 mask_opt.retrieveOption(argc,argv);
232 msknodata_opt.retrieveOption(argc,argv);
233 nodata_opt.retrieveOption(argc,argv);
235 band_opt.retrieveOption(argc,argv);
236 bstart_opt.retrieveOption(argc,argv);
237 bend_opt.retrieveOption(argc,argv);
238 balance_opt.retrieveOption(argc,argv);
239 minSize_opt.retrieveOption(argc,argv);
240 bag_opt.retrieveOption(argc,argv);
241 bagSize_opt.retrieveOption(argc,argv);
242 comb_opt.retrieveOption(argc,argv);
243 classBag_opt.retrieveOption(argc,argv);
244 prob_opt.retrieveOption(argc,argv);
245 priorimg_opt.retrieveOption(argc,argv);
246 offset_opt.retrieveOption(argc,argv);
247 scale_opt.retrieveOption(argc,argv);
248 svm_type_opt.retrieveOption(argc,argv);
249 kernel_type_opt.retrieveOption(argc,argv);
250 kernel_degree_opt.retrieveOption(argc,argv);
251 coef0_opt.retrieveOption(argc,argv);
252 nu_opt.retrieveOption(argc,argv);
253 epsilon_loss_opt.retrieveOption(argc,argv);
254 cache_opt.retrieveOption(argc,argv);
255 epsilon_tol_opt.retrieveOption(argc,argv);
256 shrinking_opt.retrieveOption(argc,argv);
257 prob_est_opt.retrieveOption(argc,argv);
258 entropy_opt.retrieveOption(argc,argv);
259 active_opt.retrieveOption(argc,argv);
260 nactive_opt.retrieveOption(argc,argv);
261 verbose_opt.retrieveOption(argc,argv);
262 random_opt.retrieveOption(argc,argv);
264 catch(
string predefinedString){
265 std::cout << predefinedString << std::endl;
270 cout <<
"Usage: pksvm -t training [-i input -o output] [-cv value]" << endl;
272 std::cout <<
"short option -h shows basic options only, use long option --help to show all options" << std::endl;
276 if(entropy_opt[0]==
"")
278 if(active_opt[0]==
"")
280 if(priorimg_opt[0]==
"")
281 priorimg_opt.clear();
284 std::map<std::string, svm::SVM_TYPE> svmMap;
286 svmMap[
"C_SVC"]=svm::C_SVC;
287 svmMap[
"nu_SVC"]=svm::nu_SVC;
288 svmMap[
"one_class"]=svm::one_class;
289 svmMap[
"epsilon_SVR"]=svm::epsilon_SVR;
290 svmMap[
"nu_SVR"]=svm::nu_SVR;
292 std::map<std::string, svm::KERNEL_TYPE> kernelMap;
294 kernelMap[
"linear"]=svm::linear;
295 kernelMap[
"polynomial"]=svm::polynomial;
296 kernelMap[
"radial"]=svm::radial;
297 kernelMap[
"sigmoid;"]=svm::sigmoid;
299 assert(training_opt.size());
301 if(verbose_opt[0]>=1){
303 std::cout <<
"input filename: " << input_opt[0] << std::endl;
305 std::cout <<
"mask filename: " << mask_opt[0] << std::endl;
306 std::cout <<
"training vector file: " << std::endl;
307 for(
int ifile=0;ifile<training_opt.size();++ifile)
308 std::cout << training_opt[ifile] << std::endl;
309 std::cout <<
"verbose: " << verbose_opt[0] << std::endl;
311 unsigned short nbag=(training_opt.size()>1)?training_opt.size():bag_opt[0];
312 if(verbose_opt[0]>=1)
313 std::cout <<
"number of bootstrap aggregations: " << nbag << std::endl;
323 bool maskIsVector=
false;
326 extentReader.open(mask_opt[0]);
328 readLayer = extentReader.getDataSource()->GetLayer(0);
329 if(!(extentReader.getExtent(ulx,uly,lrx,lry))){
330 cerr <<
"Error: could not get extent from " << mask_opt[0] << endl;
334 catch(
string errorString){
340 if(active_opt.size()){
341 prob_est_opt[0]=
true;
343 activeWriter.open(active_opt[0],ogrformat_opt[0]);
344 activeWriter.createLayer(active_opt[0],trainingReader.getProjection(),wkbPoint,NULL);
345 activeWriter.copyFields(trainingReader);
347 vector<PosValue> activePoints(nactive_opt[0]);
348 for(
int iactive=0;iactive<activePoints.size();++iactive){
349 activePoints[iactive].value=1.0;
350 activePoints[iactive].posx=0.0;
351 activePoints[iactive].posy=0.0;
354 unsigned int totalSamples=0;
355 unsigned int nactive=0;
356 vector<struct svm_model*>
svm(nbag);
357 vector<struct svm_parameter> param(nbag);
364 if(priors_opt.size()>1){
365 priors.resize(priors_opt.size());
367 for(
short iclass=0;iclass<priors_opt.size();++iclass){
368 priors[iclass]=priors_opt[iclass];
369 normPrior+=priors[iclass];
372 for(
short iclass=0;iclass<priors_opt.size();++iclass)
373 priors[iclass]/=normPrior;
378 if(bstart_opt.size()){
379 if(bend_opt.size()!=bstart_opt.size()){
380 string errorstring=
"Error: options for start and end band indexes must be provided as pairs, missing end band";
384 for(
int ipair=0;ipair<bstart_opt.size();++ipair){
385 if(bend_opt[ipair]<=bstart_opt[ipair]){
386 string errorstring=
"Error: index for end band must be smaller then start band";
389 for(
int iband=bstart_opt[ipair];iband<=bend_opt[ipair];++iband)
390 band_opt.push_back(iband);
395 cerr << error << std::endl;
400 std::sort(band_opt.begin(),band_opt.end());
402 map<string,short> classValueMap;
403 vector<std::string> nameVector;
404 if(classname_opt.size()){
405 assert(classname_opt.size()==classvalue_opt.size());
406 for(
int iclass=0;iclass<classname_opt.size();++iclass)
407 classValueMap[classname_opt[iclass]]=classvalue_opt[iclass];
412 vector< vector<double> > offset(nbag);
413 vector< vector<double> > scale(nbag);
414 map<string,Vector2d<float> > trainingMap;
415 vector< Vector2d<float> > trainingPixels;
416 vector<string> fields;
418 vector<struct svm_problem> prob(nbag);
419 vector<struct svm_node *> x_space(nbag);
421 for(
int ibag=0;ibag<nbag;++ibag){
423 if(ibag<training_opt.size()){
425 trainingPixels.clear();
426 if(verbose_opt[0]>=1)
427 std::cout <<
"reading imageVector file " << training_opt[0] << std::endl;
432 totalSamples=trainingReaderBag.readDataImageOgr(trainingMap,fields,band_opt,label_opt[0],tlayer_opt,verbose_opt[0]);
434 totalSamples=trainingReaderBag.readDataImageOgr(trainingMap,fields,0,0,label_opt[0],tlayer_opt,verbose_opt[0]);
435 if(trainingMap.size()<2){
436 string errorstring=
"Error: could not read at least two classes from training file, did you provide class labels in training sample (see option label)?";
439 trainingReaderBag.close();
442 cerr << error << std::endl;
445 catch(std::exception& e){
446 std::cerr <<
"Error: ";
447 std::cerr << e.what() << std::endl;
448 std::cerr << CPLGetLastErrorMsg() << std::endl;
452 cerr <<
"error caught" << std::endl;
459 std::cout <<
"training pixels: " << std::endl;
460 map<string,Vector2d<float> >::iterator mapit=trainingMap.begin();
461 while(mapit!=trainingMap.end()){
463 if((mapit->second).size()<minSize_opt[0]){
464 trainingMap.erase(mapit);
467 trainingPixels.push_back(mapit->second);
469 std::cout << mapit->first <<
": " << (mapit->second).size() <<
" samples" << std::endl;
473 nclass=trainingPixels.size();
474 if(classname_opt.size())
475 assert(nclass==classname_opt.size());
476 nband=trainingPixels[0][0].size()-2;
479 assert(nclass==trainingPixels.size());
480 assert(nband==trainingPixels[0][0].size()-2);
485 if(balance_opt[0]>0){
486 while(balance_opt.size()<nclass)
487 balance_opt.push_back(balance_opt.back());
491 for(
short iclass=0;iclass<nclass;++iclass){
492 if(trainingPixels[iclass].size()>balance_opt[iclass]){
493 while(trainingPixels[iclass].size()>balance_opt[iclass]){
494 int index=rand()%trainingPixels[iclass].size();
495 trainingPixels[iclass].erase(trainingPixels[iclass].begin()+index);
499 int oldsize=trainingPixels[iclass].size();
500 for(
int isample=trainingPixels[iclass].size();isample<balance_opt[iclass];++isample){
501 int index = rand()%oldsize;
502 trainingPixels[iclass].push_back(trainingPixels[iclass][index]);
505 totalSamples+=trainingPixels[iclass].size();
510 offset[ibag].resize(nband);
511 scale[ibag].resize(nband);
512 if(offset_opt.size()>1)
513 assert(offset_opt.size()==nband);
514 if(scale_opt.size()>1)
515 assert(scale_opt.size()==nband);
516 for(
int iband=0;iband<nband;++iband){
517 if(verbose_opt[0]>=1)
518 std::cout <<
"scaling for band" << iband << std::endl;
519 offset[ibag][iband]=(offset_opt.size()==1)?offset_opt[0]:offset_opt[iband];
520 scale[ibag][iband]=(scale_opt.size()==1)?scale_opt[0]:scale_opt[iband];
522 if(scale[ibag][iband]<=0){
523 float theMin=trainingPixels[0][0][iband+startBand];
524 float theMax=trainingPixels[0][0][iband+startBand];
525 for(
short iclass=0;iclass<nclass;++iclass){
526 for(
int isample=0;isample<trainingPixels[iclass].size();++isample){
527 if(theMin>trainingPixels[iclass][isample][iband+startBand])
528 theMin=trainingPixels[iclass][isample][iband+startBand];
529 if(theMax<trainingPixels[iclass][isample][iband+startBand])
530 theMax=trainingPixels[iclass][isample][iband+startBand];
533 offset[ibag][iband]=theMin+(theMax-theMin)/2.0;
534 scale[ibag][iband]=(theMax-theMin)/2.0;
535 if(verbose_opt[0]>=1){
536 std::cout <<
"Extreme image values for band " << iband <<
": [" << theMin <<
"," << theMax <<
"]" << std::endl;
537 std::cout <<
"Using offset, scale: " << offset[ibag][iband] <<
", " << scale[ibag][iband] << std::endl;
538 std::cout <<
"scaled values for band " << iband <<
": [" << (theMin-offset[ibag][iband])/scale[ibag][iband] <<
"," << (theMax-offset[ibag][iband])/scale[ibag][iband] <<
"]" << std::endl;
544 offset[ibag].resize(nband);
545 scale[ibag].resize(nband);
546 for(
int iband=0;iband<nband;++iband){
547 offset[ibag][iband]=offset[0][iband];
548 scale[ibag][iband]=scale[0][iband];
553 if(priors_opt.size()==1){
554 priors.resize(nclass);
555 for(
short iclass=0;iclass<nclass;++iclass)
556 priors[iclass]=1.0/nclass;
558 assert(priors_opt.size()==1||priors_opt.size()==nclass);
561 while(bagSize_opt.size()<nclass)
562 bagSize_opt.push_back(bagSize_opt.back());
564 if(verbose_opt[0]>=1){
565 std::cout <<
"number of bands: " << nband << std::endl;
566 std::cout <<
"number of classes: " << nclass << std::endl;
567 if(priorimg_opt.empty()){
568 std::cout <<
"priors:";
569 for(
short iclass=0;iclass<nclass;++iclass)
570 std::cout <<
" " << priors[iclass];
571 std::cout << std::endl;
574 map<string,Vector2d<float> >::iterator mapit=trainingMap.begin();
577 while(mapit!=trainingMap.end()){
578 nameVector.push_back(mapit->first);
579 if(classValueMap.size()){
581 if(classValueMap[mapit->first]>0){
582 if(cm.getClassIndex(type2string<short>(classValueMap[mapit->first]))<0){
583 cm.pushBackClassName(type2string<short>(classValueMap[mapit->first]),doSort);
587 std::cerr <<
"Error: names in classname option are not complete, please check names in training vector and make sure classvalue is > 0" << std::endl;
592 cm.pushBackClassName(mapit->first,doSort);
597 std::cerr <<
"Error: did you provide class pairs names (-c) and integer values (-r) for each class in training vector?" << std::endl;
600 if(classname_opt.empty()){
602 for(
int iclass=0;iclass<nclass;++iclass){
604 std::cout << iclass <<
" " << cm.getClass(iclass) <<
" -> " << string2type<short>(cm.getClass(iclass)) << std::endl;
605 classValueMap[cm.getClass(iclass)]=string2type<short>(cm.getClass(iclass));
617 vector< Vector2d<float> > trainingFeatures(nclass);
618 for(
short iclass=0;iclass<nclass;++iclass){
620 if(verbose_opt[0]>=1)
621 std::cout <<
"calculating features for class " << iclass << std::endl;
624 nctraining=(bagSize_opt[iclass]<100)? trainingPixels[iclass].size()/100.0*bagSize_opt[iclass] : trainingPixels[iclass].size();
627 assert(nctraining<=trainingPixels[iclass].size());
629 if(bagSize_opt[iclass]<100)
630 random_shuffle(trainingPixels[iclass].begin(),trainingPixels[iclass].end());
632 std::cout <<
"nctraining (class " << iclass <<
"): " << nctraining << std::endl;
633 trainingFeatures[iclass].resize(nctraining);
634 for(
int isample=0;isample<nctraining;++isample){
636 for(
int iband=0;iband<nband;++iband){
637 float value=trainingPixels[iclass][isample][iband+startBand];
638 trainingFeatures[iclass][isample].push_back((value-offset[ibag][iband])/scale[ibag][iband]);
641 assert(trainingFeatures[iclass].size()==nctraining);
644 unsigned int nFeatures=trainingFeatures[0][0].size();
645 if(verbose_opt[0]>=1)
646 std::cout <<
"number of features: " << nFeatures << std::endl;
647 unsigned int ntraining=0;
648 for(
short iclass=0;iclass<nclass;++iclass)
649 ntraining+=trainingFeatures[iclass].size();
650 if(verbose_opt[0]>=1)
651 std::cout <<
"training size over all classes: " << ntraining << std::endl;
653 prob[ibag].l=ntraining;
654 prob[ibag].y = Malloc(
double,prob[ibag].l);
655 prob[ibag].x = Malloc(
struct svm_node *,prob[ibag].l);
656 x_space[ibag] = Malloc(
struct svm_node,(nFeatures+1)*ntraining);
657 unsigned long int spaceIndex=0;
659 for(
short iclass=0;iclass<nclass;++iclass){
660 for(
int isample=0;isample<trainingFeatures[iclass].size();++isample){
661 prob[ibag].x[lIndex]=&(x_space[ibag][spaceIndex]);
662 for(
int ifeature=0;ifeature<nFeatures;++ifeature){
663 x_space[ibag][spaceIndex].index=ifeature+1;
664 x_space[ibag][spaceIndex].value=trainingFeatures[iclass][isample][ifeature];
667 x_space[ibag][spaceIndex++].index=-1;
668 prob[ibag].y[lIndex]=iclass;
672 assert(lIndex==prob[ibag].l);
675 param[ibag].svm_type = svmMap[svm_type_opt[0]];
676 param[ibag].kernel_type = kernelMap[kernel_type_opt[0]];
677 param[ibag].degree = kernel_degree_opt[0];
678 param[ibag].gamma = (gamma_opt[0]>0)? gamma_opt[0] : 1.0/nFeatures;
679 param[ibag].coef0 = coef0_opt[0];
680 param[ibag].nu = nu_opt[0];
681 param[ibag].cache_size = cache_opt[0];
682 param[ibag].C = ccost_opt[0];
683 param[ibag].eps = epsilon_tol_opt[0];
684 param[ibag].p = epsilon_loss_opt[0];
685 param[ibag].shrinking = (shrinking_opt[0])? 1 : 0;
686 param[ibag].probability = (prob_est_opt[0])? 1 : 0;
687 param[ibag].nr_weight = 0;
688 param[ibag].weight_label = NULL;
689 param[ibag].weight = NULL;
690 param[ibag].verbose=(verbose_opt[0]>1)?
true:
false;
693 std::cout <<
"checking parameters" << std::endl;
694 svm_check_parameter(&prob[ibag],¶m[ibag]);
696 std::cout <<
"parameters ok, training" << std::endl;
697 svm[ibag]=svm_train(&prob[ibag],¶m[ibag]);
699 std::cout <<
"SVM is now trained" << std::endl;
702 std::cout <<
"Cross validating" << std::endl;
703 double *target = Malloc(
double,prob[ibag].l);
704 svm_cross_validation(&prob[ibag],¶m[ibag],cv_opt[0],target);
705 assert(param[ibag].svm_type != EPSILON_SVR&¶m[ibag].svm_type != NU_SVR);
707 for(
int i=0;i<prob[ibag].l;i++){
708 string refClassName=nameVector[prob[ibag].y[i]];
709 string className=nameVector[target[i]];
710 if(classValueMap.size())
711 cm.incrementResult(type2string<short>(classValueMap[refClassName]),type2string<short>(classValueMap[className]),1.0/nbag);
713 cm.incrementResult(cm.getClass(prob[ibag].y[i]),cm.getClass(target[i]),1.0/nbag);
722 assert(cm.nReference());
723 cm.setFormat(cmformat_opt[0]);
724 cm.reportSE95(
false);
725 std::cout << cm << std::endl;
744 if(input_opt.empty())
747 const char* pszMessage;
748 void* pProgressArg=NULL;
749 GDALProgressFunc pfnProgress=GDALTermProgress;
752 pfnProgress(progress,pszMessage,pProgressArg);
754 bool inputIsRaster=
true;
759 if(verbose_opt[0]>=1)
760 std::cout <<
"opening image " << input_opt[0] << std::endl;
761 testImage.
open(input_opt[0]);
764 cerr << error << std::endl;
768 if(priorimg_opt.size()){
770 if(verbose_opt[0]>=1)
771 std::cout <<
"opening prior image " << priorimg_opt[0] << std::endl;
772 priorReader.
open(priorimg_opt[0]);
777 cerr << error << std::endl;
781 cerr <<
"error caught" << std::endl;
788 if(option_opt.findSubstring(
"INTERLEAVE=")==option_opt.end()){
789 string theInterleave=
"INTERLEAVE=";
791 option_opt.push_back(theInterleave);
793 vector<char> classOut(ncol);
802 if(oformat_opt.size())
803 imageType=oformat_opt[0];
805 assert(output_opt.size());
806 if(verbose_opt[0]>=1)
807 std::cout <<
"opening class image for writing output " << output_opt[0] << std::endl;
808 if(classBag_opt.size()){
809 classImageBag.
open(classBag_opt[0],ncol,nrow,nbag,GDT_Byte,imageType,option_opt);
814 classImageOut.
open(output_opt[0],ncol,nrow,1,GDT_Byte,imageType,option_opt);
818 if(colorTable_opt.size())
821 probImage.
open(prob_opt[0],ncol,nrow,nclass,GDT_Byte,imageType,option_opt);
826 if(entropy_opt.size()){
827 entropyImage.
open(entropy_opt[0],ncol,nrow,1,GDT_Byte,imageType,option_opt);
834 cerr << error << std::endl;
841 maskWriter.
open(
"/vsimem/mask.tif",ncol,nrow,1,GDT_Float32,imageType,option_opt);
845 vector<double> burnValues(1,1);
847 extentReader.close();
851 cerr << error << std::endl;
855 cerr <<
"error caught" << std::endl;
859 mask_opt.push_back(
"/vsimem/mask.tif");
864 if(verbose_opt[0]>=1)
865 std::cout <<
"opening mask image file " << mask_opt[0] << std::endl;
866 maskReader.
open(mask_opt[0]);
869 cerr << error << std::endl;
873 cerr <<
"error caught" << std::endl;
878 for(
int iline=0;iline<nrow;++iline){
879 vector<float> buffer(ncol);
880 vector<short> lineMask;
882 if(priorimg_opt.size())
883 linePrior.resize(nclass,ncol);
886 vector<float> entropy(ncol);
888 if(classBag_opt.size())
889 classBag.resize(nbag,ncol);
892 for(
int iband=0;iband<band_opt.size();++iband){
893 if(verbose_opt[0]==2)
894 std::cout <<
"reading band " << band_opt[iband] << std::endl;
895 assert(band_opt[iband]>=0);
896 assert(band_opt[iband]<testImage.
nrOfBand());
897 testImage.
readData(buffer,iline,band_opt[iband]);
898 for(
int icol=0;icol<ncol;++icol)
899 hpixel[icol].push_back(buffer[icol]);
903 for(
int iband=0;iband<nband;++iband){
904 if(verbose_opt[0]==2)
905 std::cout <<
"reading band " << iband << std::endl;
908 testImage.
readData(buffer,iline,iband);
909 for(
int icol=0;icol<ncol;++icol)
910 hpixel[icol].push_back(buffer[icol]);
914 catch(
string theError){
915 cerr <<
"Error reading " << input_opt[0] <<
": " << theError << std::endl;
919 cerr <<
"error caught" << std::endl;
922 assert(nband==hpixel[0].size());
924 std::cout <<
"used bands: " << nband << std::endl;
926 if(priorimg_opt.size()){
928 for(
short iclass=0;iclass<nclass;++iclass){
929 if(verbose_opt.size()>1)
930 std::cout <<
"Reading " << priorimg_opt[0] <<
" band " << iclass <<
" line " << iline << std::endl;
931 priorReader.
readData(linePrior[iclass],iline,iclass);
934 catch(
string theError){
935 std::cerr <<
"Error reading " << priorimg_opt[0] <<
": " << theError << std::endl;
939 cerr <<
"error caught" << std::endl;
943 double oldRowMask=-1;
945 for(
int icol=0;icol<ncol;++icol){
946 assert(hpixel[icol].size()==nband);
947 bool doClassify=
true;
953 testImage.
image2geo(icol,iline,geox,geoy);
955 if(uly>=geoy&&lry<=geoy&&ulx<=geox&&lrx>=geox){
964 testImage.
image2geo(icol,iline,geox,geoy);
965 maskReader.
geo2image(geox,geoy,colMask,rowMask);
966 colMask=
static_cast<int>(colMask);
967 rowMask=
static_cast<int>(rowMask);
968 if(rowMask>=0&&rowMask<maskReader.
nrOfRow()&&colMask>=0&&colMask<maskReader.
nrOfCol()){
969 if(static_cast<int>(rowMask)!=
static_cast<int>(oldRowMask)){
970 assert(rowMask>=0&&rowMask<maskReader.
nrOfRow());
973 maskReader.
readData(lineMask,static_cast<int>(rowMask));
975 catch(
string errorstring){
976 cerr << errorstring << endl;
980 cerr <<
"error caught" << std::endl;
986 for(
short ivalue=0;ivalue<msknodata_opt.size();++ivalue){
988 if(lineMask[colMask]==msknodata_opt[ivalue]){
989 theMask=lineMask[colMask];
1006 if(classBag_opt.size())
1007 for(
int ibag=0;ibag<nbag;++ibag)
1008 classBag[ibag][icol]=theMask;
1009 classOut[icol]=theMask;
1014 for(
int iband=0;iband<hpixel[icol].size();++iband){
1015 if(hpixel[icol][iband]){
1023 for(
short iclass=0;iclass<nclass;++iclass)
1024 probOut[iclass][icol]=0;
1026 if(classBag_opt.size())
1027 for(
int ibag=0;ibag<nbag;++ibag)
1028 classBag[ibag][icol]=nodata_opt[0];
1029 classOut[icol]=nodata_opt[0];
1032 if(verbose_opt[0]>1)
1033 std::cout <<
"begin classification " << std::endl;
1035 for(
int ibag=0;ibag<nbag;++ibag){
1036 vector<double> result(nclass);
1039 for(
int iband=0;iband<nband;++iband){
1040 x[iband].index=iband+1;
1041 x[iband].value=(hpixel[icol][iband]-offset[ibag][iband])/scale[ibag][iband];
1044 double predict_label=0;
1045 vector<float> prValues(nclass);
1047 if(!prob_est_opt[0]){
1048 predict_label = svm_predict(
svm[ibag],x);
1049 for(
short iclass=0;iclass<nclass;++iclass){
1050 if(iclass==static_cast<short>(predict_label))
1057 assert(svm_check_probability_model(
svm[ibag]));
1058 predict_label = svm_predict_probability(
svm[ibag],x,&(result[0]));
1061 if(classBag_opt.size()){
1064 classBag[ibag][icol]=0;
1067 if(priorimg_opt.size()){
1068 for(
short iclass=0;iclass<nclass;++iclass)
1069 normPrior+=linePrior[iclass][icol];
1071 for(
short iclass=0;iclass<nclass;++iclass){
1072 if(priorimg_opt.size())
1073 priors[iclass]=linePrior[iclass][icol]/normPrior;
1074 switch(comb_opt[0]){
1077 probOut[iclass][icol]+=result[iclass]*priors[iclass];
1080 probOut[iclass][icol]*=pow(static_cast<float>(priors[iclass]),static_cast<float>(1.0-nbag)/nbag)*result[iclass];
1083 if(priors[iclass]*result[iclass]>probOut[iclass][icol])
1084 probOut[iclass][icol]=priors[iclass]*result[iclass];
1087 if(classBag_opt.size()){
1093 if(result[iclass]>maxP){
1094 maxP=result[iclass];
1095 classBag[ibag][icol]=iclass;
1106 for(
short iclass=0;iclass<nclass;++iclass){
1107 if(probOut[iclass][icol]>maxBag1){
1108 maxBag1=probOut[iclass][icol];
1109 classOut[icol]=classValueMap[nameVector[iclass]];
1111 else if(probOut[iclass][icol]>maxBag2)
1112 maxBag2=probOut[iclass][icol];
1113 normBag+=probOut[iclass][icol];
1117 for(
short iclass=0;iclass<nclass;++iclass){
1118 float prv=probOut[iclass][icol];
1120 entropy[icol]-=prv*log(prv)/log(2.0);
1123 probOut[iclass][icol]=
static_cast<short>(prv+0.5);
1128 entropy[icol]/=log(static_cast<double>(nclass))/log(2.0);
1129 entropy[icol]=
static_cast<short>(100*entropy[icol]+0.5);
1130 if(active_opt.size()){
1131 if(entropy[icol]>activePoints.back().value){
1132 activePoints.back().value=entropy[icol];
1133 activePoints.back().posx=icol;
1134 activePoints.back().posy=iline;
1137 std::cout << activePoints.back().posx <<
" " << activePoints.back().posy <<
" " << activePoints.back().value << std::endl;
1142 if(classBag_opt.size())
1143 for(
int ibag=0;ibag<nbag;++ibag)
1144 classImageBag.
writeData(classBag[ibag],iline,ibag);
1145 if(prob_opt.size()){
1146 for(
short iclass=0;iclass<nclass;++iclass)
1147 probImage.
writeData(probOut[iclass],iline,iclass);
1149 if(entropy_opt.size()){
1152 classImageOut.
writeData(classOut,iline);
1153 if(!verbose_opt[0]){
1154 progress=
static_cast<float>(iline+1.0)/classImageOut.
nrOfRow();
1155 pfnProgress(progress,pszMessage,pProgressArg);
1159 if(active_opt.size()){
1160 for(
int iactive=0;iactive<activePoints.size();++iactive){
1161 std::map<string,double> pointMap;
1162 for(
int iband=0;iband<testImage.
nrOfBand();++iband){
1164 testImage.
readData(value,static_cast<int>(activePoints[iactive].posx),static_cast<int>(activePoints[iactive].posy),iband);
1167 pointMap[fs.str()]=value;
1169 pointMap[label_opt[0]]=0;
1171 testImage.
image2geo(activePoints[iactive].posx,activePoints[iactive].posy,x,y);
1172 std::string fieldname=
"id";
1173 activeWriter.addPoint(x,y,pointMap,fieldname,++nactive);
1180 if(priorimg_opt.size())
1181 priorReader.
close();
1184 if(entropy_opt.size())
1185 entropyImage.
close();
1186 if(classBag_opt.size())
1187 classImageBag.
close();
1188 classImageOut.
close();
1191 if(active_opt.size())
1192 activeWriter.close();
1194 catch(
string errorString){
1195 std::cerr <<
"Error: errorString" << std::endl;
1198 for(
int ibag=0;ibag<nbag;++ibag){
1200 svm_destroy_param(¶m[ibag]);
1203 free(x_space[ibag]);
1204 svm_free_and_destroy_model(&(
svm[ibag]));
throw this class when syntax error in command line option
void rasterizeOgr(ImgReaderOgr &ogrReader, const std::vector< double > &burnValues, const std::vector< std::string > &controlOptions=std::vector< std::string >(), const std::vector< std::string > &layernames=std::vector< std::string >())
Rasterize an OGR vector dataset using the gdal algorithm "GDALRasterizeLayers".
bool geo2image(double x, double y, double &i, double &j) const
Convert georeferenced coordinates (x and y) to image coordinates (column and row) ...
int nrOfBand(void) const
Get the number of bands of this dataset.
void close(void)
Set the memory (in MB) to cache a number of rows in memory.
int nrOfRow(void) const
Get the number of rows of this dataset.
void setColorTable(const std::string &filename, int band=0)
Set the color table using an (ASCII) file with 5 columns (value R G B alpha)
bool image2geo(double i, double j, double &x, double &y) const
Convert image coordinates (column and row) to georeferenced coordinates (x and y) ...
void readData(T &value, int col, int row, int band=0)
Read a single pixel cell value at a specific column and row for a specific band (all indices start co...
std::string getImageType() const
Get the image type (implemented as the driver description)
bool writeData(T &value, int col, int row, int band=0)
Write a single pixel cell value at a specific column and row for a specific band (all indices start c...
void copyGeoTransform(const ImgRasterGdal &imgSrc)
Copy geotransform information from another georeferenced image.
void open(const std::string &filename, const ImgReaderGdal &imgSrc, const std::vector< std::string > &options=std::vector< std::string >())
Open an image for writing, copying image attributes from a source image. Image is directly written to...
CPLErr setProjection(const std::string &projection)
Set the projection for this dataset in well known text (wkt) format.
std::string getInterleave() const
Get the band coding (interleave)
void close(void)
Close the image.
std::string getProjection(void) const
Get the projection string (deprecated, use getProjectionRef instead)
CPLErr GDALSetNoDataValue(double noDataValue, int band=0)
Set the GDAL (internal) no data value for this data set. Only a single no data value per band is supp...
void open(const std::string &filename, const GDALAccess &readMode=GA_ReadOnly)
Open an image.
int nrOfCol(void) const
Get the number of columns of this dataset.