pktools  2.6.7
Processing Kernel for geospatial data
pksvm.cc
1 /**********************************************************************
2 pksvm.cc: classify raster image using Support Vector Machine
3 Copyright (C) 2008-2017 Pieter Kempeneers
4 
5 This file is part of pktools
6 
7 pktools is free software: you can redistribute it and/or modify
8 it under the terms of the GNU General Public License as published by
9 the Free Software Foundation, either version 3 of the License, or
10 (at your option) any later version.
11 
12 pktools is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
16 
17 You should have received a copy of the GNU General Public License
18 along 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 /******************************************************************************/
118 namespace 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 
125 using namespace std;
126 
127 int 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
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.
STL namespace.
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...
Definition: ImgReaderGdal.h:95
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...
Definition: ImgWriterGdal.h:96
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...
Definition: svm.h:12
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.
Definition: ImgRasterGdal.h:98