pktools 2.6.7
Processing Kernel for geospatial data
pksvm.cc
1/**********************************************************************
2pksvm.cc: classify raster image using Support Vector Machine
3Copyright (C) 2008-2017 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 <map>
23#include <algorithm>
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"
32
33#ifdef HAVE_CONFIG_H
34#include <config.h>
35#endif
36/******************************************************************************/
118namespace svm{
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};
121}
122
123#define Malloc(type,n) (type *)malloc((n)*sizeof(type))
124
125using namespace std;
126
127int main(int argc, char *argv[])
128{
129 vector<double> priors;
130
131 //--------------------------- command line options ------------------------------------
132 Optionpk<string> input_opt("i", "input", "input image");
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)");
134 Optionpk<string> tlayer_opt("tln", "tln", "Training layer name(s)");
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)");
140 Optionpk<unsigned short> bstart_opt("sband", "startband", "Start band sequence number");
141 Optionpk<unsigned short> bend_opt("eband", "endband", "End band sequence number");
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);
146 Optionpk<unsigned short> cv_opt("cv", "cv", "N-fold cross validation mode",0);
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");
150 Optionpk<unsigned short> kernel_degree_opt("kd", "kd", "Degree in kernel function",3);
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);
160 // Optionpk<bool> weight_opt("wi", "wi", "Set the parameter C of class i to weight*C, for C_SVC",true);
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);
162 Optionpk<unsigned short> bag_opt("bag", "bag", "Number of bootstrap aggregations", 1);
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);
167 Optionpk<unsigned short> nodata_opt("nodata", "nodata", "Nodata value to put where image is masked as nodata", 0);
168 Optionpk<string> output_opt("o", "output", "Output classification image");
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)");
172 Optionpk<string> prob_opt("prob", "prob", "Probability image.");
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");
176 Optionpk<unsigned int> nactive_opt("na", "nactive", "Number of active training points",1);
177 Optionpk<string> classname_opt("c", "class", "List of class names.");
178 Optionpk<short> classvalue_opt("r", "reclass", "List of class values (use same order as in class opt).");
179 Optionpk<short> verbose_opt("v", "verbose", "Verbose level",0,2);
180
181 oformat_opt.setHide(1);
182 option_opt.setHide(1);
183 band_opt.setHide(1);
184 bstart_opt.setHide(1);
185 bend_opt.setHide(1);
186 balance_opt.setHide(1);
187 minSize_opt.setHide(1);
188 bag_opt.setHide(1);
189 bagSize_opt.setHide(1);
190 comb_opt.setHide(1);
191 classBag_opt.setHide(1);
192 prob_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);
200 nu_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);
210
211 verbose_opt.setHide(2);
212
213 bool doProcess;//stop process when program was invoked with help option (-h --help)
214 try{
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);
234 // Advanced options
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);
263 }
264 catch(string predefinedString){
265 std::cout << predefinedString << std::endl;
266 exit(0);
267 }
268 if(!doProcess){
269 cout << endl;
270 cout << "Usage: pksvm -t training [-i input -o output] [-cv value]" << endl;
271 cout << endl;
272 std::cout << "short option -h shows basic options only, use long option --help to show all options" << std::endl;
273 exit(0);//help was invoked, stop processing
274 }
275
276 if(entropy_opt[0]=="")
277 entropy_opt.clear();
278 if(active_opt[0]=="")
279 active_opt.clear();
280 if(priorimg_opt[0]=="")
281 priorimg_opt.clear();
282
283
284 std::map<std::string, svm::SVM_TYPE> svmMap;
285
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;
291
292 std::map<std::string, svm::KERNEL_TYPE> kernelMap;
293
294 kernelMap["linear"]=svm::linear;
295 kernelMap["polynomial"]=svm::polynomial;
296 kernelMap["radial"]=svm::radial;
297 kernelMap["sigmoid;"]=svm::sigmoid;
298
299 assert(training_opt.size());
300
301 if(verbose_opt[0]>=1){
302 if(input_opt.size())
303 std::cout << "input filename: " << input_opt[0] << std::endl;
304 if(mask_opt.size())
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;
310 }
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;
314
315 ImgReaderOgr extentReader;
316 OGRLayer *readLayer;
317
318 double ulx=0;
319 double uly=0;
320 double lrx=0;
321 double lry=0;
322
323 bool maskIsVector=false;
324 if(mask_opt.size()){
325 try{
326 extentReader.open(mask_opt[0]);
327 maskIsVector=true;
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;
331 exit(1);
332 }
333 }
334 catch(string errorString){
335 maskIsVector=false;
336 }
337 }
338
339 ImgWriterOgr activeWriter;
340 if(active_opt.size()){
341 prob_est_opt[0]=true;
342 ImgReaderOgr trainingReader(training_opt[0]);
343 activeWriter.open(active_opt[0],ogrformat_opt[0]);
344 activeWriter.createLayer(active_opt[0],trainingReader.getProjection(),wkbPoint,NULL);
345 activeWriter.copyFields(trainingReader);
346 }
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;
352 }
353
354 unsigned int totalSamples=0;
355 unsigned int nactive=0;
356 vector<struct svm_model*> svm(nbag);
357 vector<struct svm_parameter> param(nbag);
358
359 short nclass=0;
360 int nband=0;
361 int startBand=2;//first two bands represent X and Y pos
362
363 //normalize priors from command line
364 if(priors_opt.size()>1){//priors from argument list
365 priors.resize(priors_opt.size());
366 double normPrior=0;
367 for(short iclass=0;iclass<priors_opt.size();++iclass){
368 priors[iclass]=priors_opt[iclass];
369 normPrior+=priors[iclass];
370 }
371 //normalize
372 for(short iclass=0;iclass<priors_opt.size();++iclass)
373 priors[iclass]/=normPrior;
374 }
375
376 //convert start and end band options to vector of band indexes
377 try{
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";
381 throw(errorstring);
382 }
383 band_opt.clear();
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";
387 throw(errorstring);
388 }
389 for(int iband=bstart_opt[ipair];iband<=bend_opt[ipair];++iband)
390 band_opt.push_back(iband);
391 }
392 }
393 }
394 catch(string error){
395 cerr << error << std::endl;
396 exit(1);
397 }
398 //sort bands
399 if(band_opt.size())
400 std::sort(band_opt.begin(),band_opt.end());
401
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];
408 }
409
410 //----------------------------------- Training -------------------------------
412 vector< vector<double> > offset(nbag);
413 vector< vector<double> > scale(nbag);
414 map<string,Vector2d<float> > trainingMap;
415 vector< Vector2d<float> > trainingPixels;//[class][sample][band]
416 vector<string> fields;
417
418 vector<struct svm_problem> prob(nbag);
419 vector<struct svm_node *> x_space(nbag);
420
421 for(int ibag=0;ibag<nbag;++ibag){
422 //organize training data
423 if(ibag<training_opt.size()){//if bag contains new training pixels
424 trainingMap.clear();
425 trainingPixels.clear();
426 if(verbose_opt[0]>=1)
427 std::cout << "reading imageVector file " << training_opt[0] << std::endl;
428 try{
429 ImgReaderOgr trainingReaderBag(training_opt[ibag]);
430 if(band_opt.size())
431 //todo: when tlayer_opt is provided, readDataImageOgr does not read any layer
432 totalSamples=trainingReaderBag.readDataImageOgr(trainingMap,fields,band_opt,label_opt[0],tlayer_opt,verbose_opt[0]);
433 else
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)?";
437 throw(errorstring);
438 }
439 trainingReaderBag.close();
440 }
441 catch(string error){
442 cerr << error << std::endl;
443 exit(1);
444 }
445 catch(std::exception& e){
446 std::cerr << "Error: ";
447 std::cerr << e.what() << std::endl;
448 std::cerr << CPLGetLastErrorMsg() << std::endl;
449 exit(1);
450 }
451 catch(...){
452 cerr << "error caught" << std::endl;
453 exit(1);
454 }
455
456 //convert map to vector
457 // short iclass=0;
458 if(verbose_opt[0]>1)
459 std::cout << "training pixels: " << std::endl;
460 map<string,Vector2d<float> >::iterator mapit=trainingMap.begin();
461 while(mapit!=trainingMap.end()){
462 //delete small classes
463 if((mapit->second).size()<minSize_opt[0]){
464 trainingMap.erase(mapit);
465 continue;
466 }
467 trainingPixels.push_back(mapit->second);
468 if(verbose_opt[0]>1)
469 std::cout << mapit->first << ": " << (mapit->second).size() << " samples" << std::endl;
470 ++mapit;
471 }
472 if(!ibag){
473 nclass=trainingPixels.size();
474 if(classname_opt.size())
475 assert(nclass==classname_opt.size());
476 nband=trainingPixels[0][0].size()-2;//X and Y//trainingPixels[0][0].size();
477 }
478 else{
479 assert(nclass==trainingPixels.size());
480 assert(nband==trainingPixels[0][0].size()-2);
481 }
482
483 //do not remove outliers here: could easily be obtained through ogr2ogr -where 'B2<110' output.shp input.shp
484 //balance training data
485 if(balance_opt[0]>0){
486 while(balance_opt.size()<nclass)
487 balance_opt.push_back(balance_opt.back());
488 if(random_opt[0])
489 srand(time(NULL));
490 totalSamples=0;
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);
496 }
497 }
498 else{
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]);
503 }
504 }
505 totalSamples+=trainingPixels[iclass].size();
506 }
507 }
508
509 //set scale and offset
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];
521 //search for min and maximum
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];
531 }
532 }
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;
539 }
540 }
541 }
542 }
543 else{//use same offset and scale
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];
549 }
550 }
551
552 if(!ibag){
553 if(priors_opt.size()==1){//default: equal priors for each class
554 priors.resize(nclass);
555 for(short iclass=0;iclass<nclass;++iclass)
556 priors[iclass]=1.0/nclass;
557 }
558 assert(priors_opt.size()==1||priors_opt.size()==nclass);
559
560 //set bagsize for each class if not done already via command line
561 while(bagSize_opt.size()<nclass)
562 bagSize_opt.push_back(bagSize_opt.back());
563
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;
572 }
573 }
574 map<string,Vector2d<float> >::iterator mapit=trainingMap.begin();
575 bool doSort=true;
576 try{
577 while(mapit!=trainingMap.end()){
578 nameVector.push_back(mapit->first);
579 if(classValueMap.size()){
580 //check if name in training is covered by classname_opt (values can not be 0)
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);
584 }
585 }
586 else{
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;
588 exit(1);
589 }
590 }
591 else
592 cm.pushBackClassName(mapit->first,doSort);
593 ++mapit;
594 }
595 }
596 catch(BadConversion conversionString){
597 std::cerr << "Error: did you provide class pairs names (-c) and integer values (-r) for each class in training vector?" << std::endl;
598 exit(1);
599 }
600 if(classname_opt.empty()){
601 //std::cerr << "Warning: no class name and value pair provided for all " << nclass << " classes, using string2type<int> instead!" << std::endl;
602 for(int iclass=0;iclass<nclass;++iclass){
603 if(verbose_opt[0])
604 std::cout << iclass << " " << cm.getClass(iclass) << " -> " << string2type<short>(cm.getClass(iclass)) << std::endl;
605 classValueMap[cm.getClass(iclass)]=string2type<short>(cm.getClass(iclass));
606 }
607 }
608
609 // if(priors_opt.size()==nameVector.size()){
610 // std::cerr << "Warning: please check if priors are provided in correct order!!!" << std::endl;
611 // for(int iclass=0;iclass<nameVector.size();++iclass)
612 // std::cerr << nameVector[iclass] << " " << priors_opt[iclass] << std::endl;
613 // }
614 }//if(!ibag)
615
616 //Calculate features of training set
617 vector< Vector2d<float> > trainingFeatures(nclass);
618 for(short iclass=0;iclass<nclass;++iclass){
619 int nctraining=0;
620 if(verbose_opt[0]>=1)
621 std::cout << "calculating features for class " << iclass << std::endl;
622 if(random_opt[0])
623 srand(time(NULL));
624 nctraining=(bagSize_opt[iclass]<100)? trainingPixels[iclass].size()/100.0*bagSize_opt[iclass] : trainingPixels[iclass].size();//bagSize_opt[iclass] given in % of training size
625 if(nctraining<=0)
626 nctraining=1;
627 assert(nctraining<=trainingPixels[iclass].size());
628 int index=0;
629 if(bagSize_opt[iclass]<100)
630 random_shuffle(trainingPixels[iclass].begin(),trainingPixels[iclass].end());
631 if(verbose_opt[0]>1)
632 std::cout << "nctraining (class " << iclass << "): " << nctraining << std::endl;
633 trainingFeatures[iclass].resize(nctraining);
634 for(int isample=0;isample<nctraining;++isample){
635 //scale pixel values according to scale and offset!!!
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]);
639 }
640 }
641 assert(trainingFeatures[iclass].size()==nctraining);
642 }
643
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;
652
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;
658 int lIndex=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];
665 ++spaceIndex;
666 }
667 x_space[ibag][spaceIndex++].index=-1;
668 prob[ibag].y[lIndex]=iclass;
669 ++lIndex;
670 }
671 }
672 assert(lIndex==prob[ibag].l);
673
674 //set SVM parameters through command line options
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;//not used: I use priors and balancing
688 param[ibag].weight_label = NULL;
689 param[ibag].weight = NULL;
690 param[ibag].verbose=(verbose_opt[0]>1)? true:false;
691
692 if(verbose_opt[0]>1)
693 std::cout << "checking parameters" << std::endl;
694 svm_check_parameter(&prob[ibag],&param[ibag]);
695 if(verbose_opt[0])
696 std::cout << "parameters ok, training" << std::endl;
697 svm[ibag]=svm_train(&prob[ibag],&param[ibag]);
698 if(verbose_opt[0]>1)
699 std::cout << "SVM is now trained" << std::endl;
700 if(cv_opt[0]>1){
701 if(verbose_opt[0]>1)
702 std::cout << "Cross validating" << std::endl;
703 double *target = Malloc(double,prob[ibag].l);
704 svm_cross_validation(&prob[ibag],&param[ibag],cv_opt[0],target);
705 assert(param[ibag].svm_type != EPSILON_SVR&&param[ibag].svm_type != NU_SVR);//only for regression
706
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);
712 else
713 cm.incrementResult(cm.getClass(prob[ibag].y[i]),cm.getClass(target[i]),1.0/nbag);
714 }
715 free(target);
716 }
717 // *NOTE* Because svm_model contains pointers to svm_problem, you can
718 // not free the memory used by svm_problem if you are still using the
719 // svm_model produced by svm_train().
720 }//for ibag
721 if(cv_opt[0]>1){
722 assert(cm.nReference());
723 cm.setFormat(cmformat_opt[0]);
724 cm.reportSE95(false);
725 std::cout << cm << std::endl;
726 // cout << "class #samples userAcc prodAcc" << endl;
727 // double se95_ua=0;
728 // double se95_pa=0;
729 // double se95_oa=0;
730 // double dua=0;
731 // double dpa=0;
732 // double doa=0;
733 // for(short iclass=0;iclass<cm.nClasses();++iclass){
734 // dua=cm.ua(cm.getClass(iclass),&se95_ua);
735 // dpa=cm.pa(cm.getClass(iclass),&se95_pa);
736 // cout << cm.getClass(iclass) << " " << cm.nReference(cm.getClass(iclass)) << " " << dua << " (" << se95_ua << ")" << " " << dpa << " (" << se95_pa << ")" << endl;
737 // }
738 // std::cout << "Kappa: " << cm.kappa() << std::endl;
739 // doa=cm.oa(&se95_oa);
740 // std::cout << "Overall Accuracy: " << 100*doa << " (" << 100*se95_oa << ")" << std::endl;
741 }
742
743 //--------------------------------- end of training -----------------------------------
744 if(input_opt.empty())
745 exit(0);
746
747 const char* pszMessage;
748 void* pProgressArg=NULL;
749 GDALProgressFunc pfnProgress=GDALTermProgress;
750 float progress=0;
751 if(!verbose_opt[0])
752 pfnProgress(progress,pszMessage,pProgressArg);
753 //-------------------------------- open image file ------------------------------------
754 bool inputIsRaster=true;
755 ImgReaderOgr imgReaderOgr;
756 if(inputIsRaster){
757 ImgReaderGdal testImage;
758 try{
759 if(verbose_opt[0]>=1)
760 std::cout << "opening image " << input_opt[0] << std::endl;
761 testImage.open(input_opt[0]);
762 }
763 catch(string error){
764 cerr << error << std::endl;
765 exit(2);
766 }
767 ImgReaderGdal priorReader;
768 if(priorimg_opt.size()){
769 try{
770 if(verbose_opt[0]>=1)
771 std::cout << "opening prior image " << priorimg_opt[0] << std::endl;
772 priorReader.open(priorimg_opt[0]);
773 assert(priorReader.nrOfCol()==testImage.nrOfCol());
774 assert(priorReader.nrOfRow()==testImage.nrOfRow());
775 }
776 catch(string error){
777 cerr << error << std::endl;
778 exit(2);
779 }
780 catch(...){
781 cerr << "error caught" << std::endl;
782 exit(1);
783 }
784 }
785
786 int nrow=testImage.nrOfRow();
787 int ncol=testImage.nrOfCol();
788 if(option_opt.findSubstring("INTERLEAVE=")==option_opt.end()){
789 string theInterleave="INTERLEAVE=";
790 theInterleave+=testImage.getInterleave();
791 option_opt.push_back(theInterleave);
792 }
793 vector<char> classOut(ncol);//classified line for writing to image file
794
795 // assert(nband==testImage.nrOfBand());
796 ImgWriterGdal classImageBag;
797 ImgWriterGdal classImageOut;
798 ImgWriterGdal probImage;
799 ImgWriterGdal entropyImage;
800
801 string imageType=testImage.getImageType();
802 if(oformat_opt.size())//default
803 imageType=oformat_opt[0];
804 try{
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);
810 classImageBag.GDALSetNoDataValue(nodata_opt[0]);
811 classImageBag.copyGeoTransform(testImage);
812 classImageBag.setProjection(testImage.getProjection());
813 }
814 classImageOut.open(output_opt[0],ncol,nrow,1,GDT_Byte,imageType,option_opt);
815 classImageOut.GDALSetNoDataValue(nodata_opt[0]);
816 classImageOut.copyGeoTransform(testImage);
817 classImageOut.setProjection(testImage.getProjection());
818 if(colorTable_opt.size())
819 classImageOut.setColorTable(colorTable_opt[0],0);
820 if(prob_opt.size()){
821 probImage.open(prob_opt[0],ncol,nrow,nclass,GDT_Byte,imageType,option_opt);
822 probImage.GDALSetNoDataValue(nodata_opt[0]);
823 probImage.copyGeoTransform(testImage);
824 probImage.setProjection(testImage.getProjection());
825 }
826 if(entropy_opt.size()){
827 entropyImage.open(entropy_opt[0],ncol,nrow,1,GDT_Byte,imageType,option_opt);
828 entropyImage.GDALSetNoDataValue(nodata_opt[0]);
829 entropyImage.copyGeoTransform(testImage);
830 entropyImage.setProjection(testImage.getProjection());
831 }
832 }
833 catch(string error){
834 cerr << error << std::endl;
835 }
836
837 ImgWriterGdal maskWriter;
838
839 if(maskIsVector){
840 try{
841 maskWriter.open("/vsimem/mask.tif",ncol,nrow,1,GDT_Float32,imageType,option_opt);
842 maskWriter.GDALSetNoDataValue(nodata_opt[0]);
843 maskWriter.copyGeoTransform(testImage);
844 maskWriter.setProjection(testImage.getProjection());
845 vector<double> burnValues(1,1);//burn value is 1 (single band)
846 maskWriter.rasterizeOgr(extentReader,burnValues);
847 extentReader.close();
848 maskWriter.close();
849 }
850 catch(string error){
851 cerr << error << std::endl;
852 exit(2);
853 }
854 catch(...){
855 cerr << "error caught" << std::endl;
856 exit(1);
857 }
858 mask_opt.clear();
859 mask_opt.push_back("/vsimem/mask.tif");
860 }
861 ImgReaderGdal maskReader;
862 if(mask_opt.size()){
863 try{
864 if(verbose_opt[0]>=1)
865 std::cout << "opening mask image file " << mask_opt[0] << std::endl;
866 maskReader.open(mask_opt[0]);
867 }
868 catch(string error){
869 cerr << error << std::endl;
870 exit(2);
871 }
872 catch(...){
873 cerr << "error caught" << std::endl;
874 exit(1);
875 }
876 }
877
878 for(int iline=0;iline<nrow;++iline){
879 vector<float> buffer(ncol);
880 vector<short> lineMask;
881 Vector2d<float> linePrior;
882 if(priorimg_opt.size())
883 linePrior.resize(nclass,ncol);//prior prob for each class
884 Vector2d<float> hpixel(ncol);
885 Vector2d<float> probOut(nclass,ncol);//posterior prob for each (internal) class
886 vector<float> entropy(ncol);
887 Vector2d<char> classBag;//classified line for writing to image file
888 if(classBag_opt.size())
889 classBag.resize(nbag,ncol);
890 try{
891 if(band_opt.size()){
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]);
900 }
901 }
902 else{
903 for(int iband=0;iband<nband;++iband){
904 if(verbose_opt[0]==2)
905 std::cout << "reading band " << iband << std::endl;
906 assert(iband>=0);
907 assert(iband<testImage.nrOfBand());
908 testImage.readData(buffer,iline,iband);
909 for(int icol=0;icol<ncol;++icol)
910 hpixel[icol].push_back(buffer[icol]);
911 }
912 }
913 }
914 catch(string theError){
915 cerr << "Error reading " << input_opt[0] << ": " << theError << std::endl;
916 exit(3);
917 }
918 catch(...){
919 cerr << "error caught" << std::endl;
920 exit(3);
921 }
922 assert(nband==hpixel[0].size());
923 if(verbose_opt[0]>1)
924 std::cout << "used bands: " << nband << std::endl;
925 //read prior
926 if(priorimg_opt.size()){
927 try{
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);
932 }
933 }
934 catch(string theError){
935 std::cerr << "Error reading " << priorimg_opt[0] << ": " << theError << std::endl;
936 exit(3);
937 }
938 catch(...){
939 cerr << "error caught" << std::endl;
940 exit(3);
941 }
942 }
943 double oldRowMask=-1;//keep track of row mask to optimize number of line readings
944 //process per pixel
945 for(int icol=0;icol<ncol;++icol){
946 assert(hpixel[icol].size()==nband);
947 bool doClassify=true;
948 bool masked=false;
949 double geox=0;
950 double geoy=0;
951 if(maskIsVector){
952 doClassify=false;
953 testImage.image2geo(icol,iline,geox,geoy);
954 //check enveloppe first
955 if(uly>=geoy&&lry<=geoy&&ulx<=geox&&lrx>=geox){
956 doClassify=true;
957 }
958 }
959 if(mask_opt.size()){
960 //read mask
961 double colMask=0;
962 double rowMask=0;
963
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());
971 try{
972 // maskReader.readData(lineMask[imask],GDT_Int32,static_cast<int>(rowMask));
973 maskReader.readData(lineMask,static_cast<int>(rowMask));
974 }
975 catch(string errorstring){
976 cerr << errorstring << endl;
977 exit(1);
978 }
979 catch(...){
980 cerr << "error caught" << std::endl;
981 exit(3);
982 }
983 oldRowMask=rowMask;
984 }
985 short theMask=0;
986 for(short ivalue=0;ivalue<msknodata_opt.size();++ivalue){
987 // if(msknodata_opt[ivalue]>=0){//values set in msknodata_opt are invalid
988 if(lineMask[colMask]==msknodata_opt[ivalue]){
989 theMask=lineMask[colMask];
990 masked=true;
991 break;
992 }
993 // }
994 // else{//only values set in msknodata_opt are valid
995 // if(lineMask[colMask]!=-msknodata_opt[ivalue]){
996 // theMask=lineMask[colMask];
997 // masked=true;
998 // }
999 // else{
1000 // masked=false;
1001 // break;
1002 // }
1003 // }
1004 }
1005 if(masked){
1006 if(classBag_opt.size())
1007 for(int ibag=0;ibag<nbag;++ibag)
1008 classBag[ibag][icol]=theMask;
1009 classOut[icol]=theMask;
1010 continue;
1011 }
1012 }
1013 bool valid=false;
1014 for(int iband=0;iband<hpixel[icol].size();++iband){
1015 if(hpixel[icol][iband]){
1016 valid=true;
1017 break;
1018 }
1019 }
1020 if(!valid)
1021 doClassify=false;
1022 }
1023 for(short iclass=0;iclass<nclass;++iclass)
1024 probOut[iclass][icol]=0;
1025 if(!doClassify){
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];
1030 continue;//next column
1031 }
1032 if(verbose_opt[0]>1)
1033 std::cout << "begin classification " << std::endl;
1034 //----------------------------------- classification -------------------
1035 for(int ibag=0;ibag<nbag;++ibag){
1036 vector<double> result(nclass);
1037 struct svm_node *x;
1038 x = (struct svm_node *) malloc((nband+1)*sizeof(struct svm_node));
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];
1042 }
1043 x[nband].index=-1;//to end svm feature vector
1044 double predict_label=0;
1045 vector<float> prValues(nclass);
1046 float maxP=0;
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))
1051 result[iclass]=1;
1052 else
1053 result[iclass]=0;
1054 }
1055 }
1056 else{
1057 assert(svm_check_probability_model(svm[ibag]));
1058 predict_label = svm_predict_probability(svm[ibag],x,&(result[0]));
1059 }
1060 //calculate posterior prob of bag
1061 if(classBag_opt.size()){
1062 //search for max prob within bag
1063 maxP=0;
1064 classBag[ibag][icol]=0;
1065 }
1066 double normPrior=0;
1067 if(priorimg_opt.size()){
1068 for(short iclass=0;iclass<nclass;++iclass)
1069 normPrior+=linePrior[iclass][icol];
1070 }
1071 for(short iclass=0;iclass<nclass;++iclass){
1072 if(priorimg_opt.size())
1073 priors[iclass]=linePrior[iclass][icol]/normPrior;//todo: check if correct for all cases... (automatic classValueMap and manual input for names and values)
1074 switch(comb_opt[0]){
1075 default:
1076 case(0)://sum rule
1077 probOut[iclass][icol]+=result[iclass]*priors[iclass];//add probabilities for each bag
1078 break;
1079 case(1)://product rule
1080 probOut[iclass][icol]*=pow(static_cast<float>(priors[iclass]),static_cast<float>(1.0-nbag)/nbag)*result[iclass];//multiply probabilities for each bag
1081 break;
1082 case(2)://max rule
1083 if(priors[iclass]*result[iclass]>probOut[iclass][icol])
1084 probOut[iclass][icol]=priors[iclass]*result[iclass];
1085 break;
1086 }
1087 if(classBag_opt.size()){
1088 //search for max prob within bag
1089 // if(prValues[iclass]>maxP){
1090 // maxP=prValues[iclass];
1091 // classBag[ibag][icol]=iclass;
1092 // }
1093 if(result[iclass]>maxP){
1094 maxP=result[iclass];
1095 classBag[ibag][icol]=iclass;
1096 }
1097 }
1098 }
1099 free(x);
1100 }//ibag
1101
1102 //search for max class prob
1103 float maxBag1=0;//max probability
1104 float maxBag2=0;//second max probability
1105 float normBag=0;
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]];
1110 }
1111 else if(probOut[iclass][icol]>maxBag2)
1112 maxBag2=probOut[iclass][icol];
1113 normBag+=probOut[iclass][icol];
1114 }
1115 //normalize probOut and convert to percentage
1116 entropy[icol]=0;
1117 for(short iclass=0;iclass<nclass;++iclass){
1118 float prv=probOut[iclass][icol];
1119 prv/=normBag;
1120 entropy[icol]-=prv*log(prv)/log(2.0);
1121 prv*=100.0;
1122
1123 probOut[iclass][icol]=static_cast<short>(prv+0.5);
1124 // assert(classValueMap[nameVector[iclass]]<probOut.size());
1125 // assert(classValueMap[nameVector[iclass]]>=0);
1126 // probOut[classValueMap[nameVector[iclass]]][icol]=static_cast<short>(prv+0.5);
1127 }
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];//replace largest value (last)
1133 activePoints.back().posx=icol;
1134 activePoints.back().posy=iline;
1135 std::sort(activePoints.begin(),activePoints.end(),Decrease_PosValue());//sort in descending order (largest first, smallest last)
1136 if(verbose_opt[0])
1137 std::cout << activePoints.back().posx << " " << activePoints.back().posy << " " << activePoints.back().value << std::endl;
1138 }
1139 }
1140 }//icol
1141 //----------------------------------- write output ------------------------------------------
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);
1148 }
1149 if(entropy_opt.size()){
1150 entropyImage.writeData(entropy,iline);
1151 }
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);
1156 }
1157 }
1158 //write active learning points
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){
1163 double value;
1164 testImage.readData(value,static_cast<int>(activePoints[iactive].posx),static_cast<int>(activePoints[iactive].posy),iband);
1165 ostringstream fs;
1166 fs << "B" << iband;
1167 pointMap[fs.str()]=value;
1168 }
1169 pointMap[label_opt[0]]=0;
1170 double x, y;
1171 testImage.image2geo(activePoints[iactive].posx,activePoints[iactive].posy,x,y);
1172 std::string fieldname="id";//number of the point
1173 activeWriter.addPoint(x,y,pointMap,fieldname,++nactive);
1174 }
1175 }
1176
1177 testImage.close();
1178 if(mask_opt.size())
1179 maskReader.close();
1180 if(priorimg_opt.size())
1181 priorReader.close();
1182 if(prob_opt.size())
1183 probImage.close();
1184 if(entropy_opt.size())
1185 entropyImage.close();
1186 if(classBag_opt.size())
1187 classImageBag.close();
1188 classImageOut.close();
1189 }
1190 try{
1191 if(active_opt.size())
1192 activeWriter.close();
1193 }
1194 catch(string errorString){
1195 std::cerr << "Error: errorString" << std::endl;
1196 }
1197
1198 for(int ibag=0;ibag<nbag;++ibag){
1199 // svm_destroy_param[ibag](&param[ibag]);
1200 svm_destroy_param(&param[ibag]);
1201 free(prob[ibag].y);
1202 free(prob[ibag].x);
1203 free(x_space[ibag]);
1204 svm_free_and_destroy_model(&(svm[ibag]));
1205 }
1206 return 0;
1207}
throw this class when syntax error in command line option
Definition: Optionpk.h:45
int nrOfRow(void) const
Get the number of rows of this dataset.
bool geo2image(double x, double y, double &i, double &j) const
Convert georeferenced coordinates (x and y) to image coordinates (column and row)
std::string getProjection(void) const
Get the projection string (deprecated, use getProjectionRef instead)
int nrOfBand(void) const
Get the number of bands of this dataset.
void copyGeoTransform(const ImgRasterGdal &imgSrc)
Copy geotransform information from another georeferenced image.
int nrOfCol(void) const
Get the number of columns of this dataset.
Definition: ImgRasterGdal.h:98
std::string getInterleave() const
Get the band coding (interleave)
bool image2geo(double i, double j, double &x, double &y) const
Convert image coordinates (column and row) to georeferenced coordinates (x and y)
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...
std::string getImageType() const
Get the image type (implemented as the driver description)
CPLErr setProjection(const std::string &projection)
Set the projection for this dataset in well known text (wkt) format.
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...
Definition: ImgReaderGdal.h:95
void close(void)
Set the memory (in MB) to cache a number of rows in memory.
void open(const std::string &filename, const GDALAccess &readMode=GA_ReadOnly)
Open an image.
void close(void)
Close the 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...
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".
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 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...
Definition: ImgWriterGdal.h:96
Definition: svm.h:13