pktools 2.6.7
Processing Kernel for geospatial data
pkannogr.cc
1/**********************************************************************
2pkannogr.cc: classify vector dataset using Artificial Neural Network
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 "floatfann.h"
32#include "algorithms/myfann_cpp.h"
33
34/******************************************************************************/
104using namespace std;
105
106int main(int argc, char *argv[])
107{
108 vector<double> priors;
109
110 //--------------------------- command line options ------------------------------------
111 Optionpk<string> input_opt("i", "input", "input image");
112 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)");
113 Optionpk<string> tlayer_opt("tln", "tln", "training layer name(s)");
114 Optionpk<string> label_opt("label", "label", "identifier for class label in training vector file.","label");
115 Optionpk<unsigned int> balance_opt("bal", "balance", "balance the input data to this number of samples for each class", 0);
116 Optionpk<bool> random_opt("random", "random", "in case of balance, randomize input data", true,2);
117 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);
118 Optionpk<unsigned short> band_opt("b", "band", "band index (starting from 0, either use band option or use start to end)");
119 Optionpk<unsigned short> bstart_opt("sband", "startband", "Start band sequence number");
120 Optionpk<unsigned short> bend_opt("eband", "endband", "End band sequence number");
121 Optionpk<double> offset_opt("offset", "offset", "offset value for each spectral band input features: refl[band]=(DN[band]-offset[band])/scale[band]", 0.0);
122 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);
123 Optionpk<unsigned short> aggreg_opt("a", "aggreg", "how to combine aggregated classifiers, see also rc option (1: sum rule, 2: max rule).",1);
124 Optionpk<double> priors_opt("prior", "prior", "prior probabilities for each class (e.g., -p 0.3 -p 0.3 -p 0.2 )", 0.0);
125 Optionpk<string> priorimg_opt("pim", "priorimg", "prior probability image (multi-band img with band for each class","",2);
126 Optionpk<unsigned short> cv_opt("cv", "cv", "n-fold cross validation mode",0);
127 Optionpk<string> cmformat_opt("cmf","cmf","Format for confusion matrix (ascii or latex)","ascii");
128 Optionpk<unsigned int> nneuron_opt("nn", "nneuron", "number of neurons in hidden layers in neural network (multiple hidden layers are set by defining multiple number of neurons: -n 15 -n 1, default is one hidden layer with 5 neurons)", 5);
129 Optionpk<float> connection_opt("\0", "connection", "connection reate (default: 1.0 for a fully connected network)", 1.0);
130 Optionpk<float> learning_opt("l", "learning", "learning rate (default: 0.7)", 0.7);
131 Optionpk<float> weights_opt("w", "weights", "weights for neural network. Apply to fully connected network only, starting from first input neuron to last output neuron, including the bias neurons (last neuron in each but last layer)", 0.0);
132 Optionpk<unsigned int> maxit_opt("\0", "maxit", "number of maximum iterations (epoch) (default: 500)", 500);
133 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. Default is sum rule (0)",0);
134 Optionpk<unsigned short> bag_opt("bag", "bag", "Number of bootstrap aggregations (default is no bagging: 1)", 1);
135 Optionpk<int> bagSize_opt("bs", "bsize", "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);
136 Optionpk<string> classBag_opt("cb", "classbag", "output for each individual bootstrap aggregation (default is blank)");
137 Optionpk<string> mask_opt("m", "mask", "Only classify within specified mask (vector or raster).");
138 Optionpk<unsigned short> nodata_opt("nodata", "nodata", "nodata value to put where image is masked as nodata", 0);
139 Optionpk<string> output_opt("o", "output", "output classification image");
140 Optionpk<string> otype_opt("ot", "otype", "Data type for output image ({Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/CInt16/CInt32/CFloat32/CFloat64}). Empty string: inherit type from input image");
141 Optionpk<string> option_opt("co", "co", "Creation option for output file. Multiple options can be specified.");
142 Optionpk<string> colorTable_opt("ct", "ct", "colour table in ASCII format having 5 columns: id R G B ALFA (0: transparent, 255: solid)");
143 Optionpk<string> entropy_opt("entropy", "entropy", "entropy image (measure for uncertainty of classifier output","",2);
144 Optionpk<string> active_opt("active", "active", "ogr output for active training sample.","",2);
145 Optionpk<string> ogrformat_opt("f", "f", "Output ogr format for active training sample","SQLite");
146 Optionpk<unsigned int> nactive_opt("na", "nactive", "number of active training points",1);
147 Optionpk<string> classname_opt("c", "class", "list of class names.");
148 Optionpk<short> classvalue_opt("r", "reclass", "list of class values (use same order as in class opt).");
149 Optionpk<short> verbose_opt("v", "verbose", "set to: 0 (results only), 1 (confusion matrix), 2 (debug)",0,2);
150
151 option_opt.setHide(1);
152 band_opt.setHide(1);
153 bstart_opt.setHide(1);
154 bend_opt.setHide(1);
155 balance_opt.setHide(1);
156 minSize_opt.setHide(1);
157 bag_opt.setHide(1);
158 bagSize_opt.setHide(1);
159 comb_opt.setHide(1);
160 classBag_opt.setHide(1);
161 minSize_opt.setHide(1);
162 priorimg_opt.setHide(1);
163 minSize_opt.setHide(1);
164 offset_opt.setHide(1);
165 scale_opt.setHide(1);
166 connection_opt.setHide(1);
167 weights_opt.setHide(1);
168 maxit_opt.setHide(1);
169 learning_opt.setHide(1);
170
171 verbose_opt.setHide(2);
172
173 bool doProcess;//stop process when program was invoked with help option (-h --help)
174 try{
175 doProcess=input_opt.retrieveOption(argc,argv);
176 training_opt.retrieveOption(argc,argv);
177 tlayer_opt.retrieveOption(argc,argv);
178 label_opt.retrieveOption(argc,argv);
179 balance_opt.retrieveOption(argc,argv);
180 random_opt.retrieveOption(argc,argv);
181 minSize_opt.retrieveOption(argc,argv);
182 band_opt.retrieveOption(argc,argv);
183 bstart_opt.retrieveOption(argc,argv);
184 bend_opt.retrieveOption(argc,argv);
185 offset_opt.retrieveOption(argc,argv);
186 scale_opt.retrieveOption(argc,argv);
187 aggreg_opt.retrieveOption(argc,argv);
188 priors_opt.retrieveOption(argc,argv);
189 priorimg_opt.retrieveOption(argc,argv);
190 cv_opt.retrieveOption(argc,argv);
191 cmformat_opt.retrieveOption(argc,argv);
192 nneuron_opt.retrieveOption(argc,argv);
193 connection_opt.retrieveOption(argc,argv);
194 weights_opt.retrieveOption(argc,argv);
195 learning_opt.retrieveOption(argc,argv);
196 maxit_opt.retrieveOption(argc,argv);
197 comb_opt.retrieveOption(argc,argv);
198 bag_opt.retrieveOption(argc,argv);
199 bagSize_opt.retrieveOption(argc,argv);
200 classBag_opt.retrieveOption(argc,argv);
201 mask_opt.retrieveOption(argc,argv);
202 nodata_opt.retrieveOption(argc,argv);
203 output_opt.retrieveOption(argc,argv);
204 otype_opt.retrieveOption(argc,argv);
205 colorTable_opt.retrieveOption(argc,argv);
206 option_opt.retrieveOption(argc,argv);
207 entropy_opt.retrieveOption(argc,argv);
208 active_opt.retrieveOption(argc,argv);
209 ogrformat_opt.retrieveOption(argc,argv);
210 nactive_opt.retrieveOption(argc,argv);
211 classname_opt.retrieveOption(argc,argv);
212 classvalue_opt.retrieveOption(argc,argv);
213 verbose_opt.retrieveOption(argc,argv);
214 }
215 catch(string predefinedString){
216 std::cout << predefinedString << std::endl;
217 exit(0);
218 }
219 if(!doProcess){
220 cout << endl;
221 cout << "Usage: pkannogr -t training [-i input -o output] [-cv value]" << endl;
222 cout << endl;
223 cout << "short option -h shows basic options only, use long option --help to show all options" << endl;
224 exit(0);//help was invoked, stop processing
225 }
226
227 if(entropy_opt[0]=="")
228 entropy_opt.clear();
229 if(active_opt[0]=="")
230 active_opt.clear();
231 if(priorimg_opt[0]=="")
232 priorimg_opt.clear();
233
234 if(verbose_opt[0]>=1){
235 if(input_opt.size())
236 cout << "image filename: " << input_opt[0] << endl;
237 if(mask_opt.size())
238 cout << "mask filename: " << mask_opt[0] << endl;
239 if(training_opt.size()){
240 cout << "training vector file: " << endl;
241 for(int ifile=0;ifile<training_opt.size();++ifile)
242 cout << training_opt[ifile] << endl;
243 }
244 else
245 cerr << "no training file set!" << endl;
246 cout << "verbose: " << verbose_opt[0] << endl;
247 }
248 unsigned short nbag=(training_opt.size()>1)?training_opt.size():bag_opt[0];
249 if(verbose_opt[0]>=1)
250 cout << "number of bootstrap aggregations: " << nbag << endl;
251
252 ImgReaderOgr extentReader;
253 OGRLayer *readLayer;
254
255 double ulx=0;
256 double uly=0;
257 double lrx=0;
258 double lry=0;
259
260 bool maskIsVector=false;
261 if(mask_opt.size()){
262 try{
263 extentReader.open(mask_opt[0]);
264 maskIsVector=true;
265 readLayer = extentReader.getDataSource()->GetLayer(0);
266 if(!(extentReader.getExtent(ulx,uly,lrx,lry))){
267 cerr << "Error: could not get extent from " << mask_opt[0] << endl;
268 exit(1);
269 }
270 }
271 catch(string errorString){
272 maskIsVector=false;
273 }
274 }
275
276 ImgWriterOgr activeWriter;
277 if(active_opt.size()){
278 ImgReaderOgr trainingReader(training_opt[0]);
279 activeWriter.open(active_opt[0],ogrformat_opt[0]);
280 activeWriter.createLayer(active_opt[0],trainingReader.getProjection(),wkbPoint,NULL);
281 activeWriter.copyFields(trainingReader);
282 }
283 vector<PosValue> activePoints(nactive_opt[0]);
284 for(int iactive=0;iactive<activePoints.size();++iactive){
285 activePoints[iactive].value=1.0;
286 activePoints[iactive].posx=0.0;
287 activePoints[iactive].posy=0.0;
288 }
289
290 unsigned int totalSamples=0;
291 unsigned int nactive=0;
292 vector<FANN::neural_net> net(nbag);//the neural network
293
294 unsigned int nclass=0;
295 int nband=0;
296 int startBand=2;//first two bands represent X and Y pos
297
298 if(priors_opt.size()>1){//priors from argument list
299 priors.resize(priors_opt.size());
300 double normPrior=0;
301 for(int iclass=0;iclass<priors_opt.size();++iclass){
302 priors[iclass]=priors_opt[iclass];
303 normPrior+=priors[iclass];
304 }
305 //normalize
306 for(int iclass=0;iclass<priors_opt.size();++iclass)
307 priors[iclass]/=normPrior;
308 }
309
310 //convert start and end band options to vector of band indexes
311 try{
312 if(bstart_opt.size()){
313 if(bend_opt.size()!=bstart_opt.size()){
314 string errorstring="Error: options for start and end band indexes must be provided as pairs, missing end band";
315 throw(errorstring);
316 }
317 band_opt.clear();
318 for(int ipair=0;ipair<bstart_opt.size();++ipair){
319 if(bend_opt[ipair]<=bstart_opt[ipair]){
320 string errorstring="Error: index for end band must be smaller then start band";
321 throw(errorstring);
322 }
323 for(int iband=bstart_opt[ipair];iband<=bend_opt[ipair];++iband)
324 band_opt.push_back(iband);
325 }
326 }
327 }
328 catch(string error){
329 cerr << error << std::endl;
330 exit(1);
331 }
332 //sort bands
333 if(band_opt.size())
334 std::sort(band_opt.begin(),band_opt.end());
335
336 map<string,short> classValueMap;
337 vector<std::string> nameVector;
338 if(classname_opt.size()){
339 assert(classname_opt.size()==classvalue_opt.size());
340 for(int iclass=0;iclass<classname_opt.size();++iclass)
341 classValueMap[classname_opt[iclass]]=classvalue_opt[iclass];
342 }
343 //----------------------------------- Training -------------------------------
345 vector< vector<double> > offset(nbag);
346 vector< vector<double> > scale(nbag);
347 map<string,Vector2d<float> > trainingMap;
348 vector< Vector2d<float> > trainingPixels;//[class][sample][band]
349 vector<string> fields;
350 for(int ibag=0;ibag<nbag;++ibag){
351 //organize training data
352 if(ibag<training_opt.size()){//if bag contains new training pixels
353 trainingMap.clear();
354 trainingPixels.clear();
355 if(verbose_opt[0]>=1)
356 cout << "reading imageVector file " << training_opt[0] << endl;
357 try{
358 ImgReaderOgr trainingReaderBag(training_opt[ibag]);
359 if(band_opt.size())
360 totalSamples=trainingReaderBag.readDataImageOgr(trainingMap,fields,band_opt,label_opt[0],tlayer_opt);
361 else
362 totalSamples=trainingReaderBag.readDataImageOgr(trainingMap,fields,0,0,label_opt[0],tlayer_opt);
363 if(trainingMap.size()<2){
364 string errorstring="Error: could not read at least two classes from training file, did you provide class labels in training sample (see option label)?";
365 throw(errorstring);
366 }
367 trainingReaderBag.close();
368 }
369 catch(string error){
370 cerr << error << std::endl;
371 exit(1);
372 }
373 catch(...){
374 cerr << "error caught" << std::endl;
375 exit(1);
376 }
377 //delete class 0 ?
378 // if(verbose_opt[0]>=1)
379 // std::cout << "erasing class 0 from training set (" << trainingMap[0].size() << " from " << totalSamples << ") samples" << std::endl;
380 // totalSamples-=trainingMap[0].size();
381 // trainingMap.erase(0);
382 //convert map to vector
383 if(verbose_opt[0]>1)
384 std::cout << "training pixels: " << std::endl;
385 map<string,Vector2d<float> >::iterator mapit=trainingMap.begin();
386 while(mapit!=trainingMap.end()){
387 //delete small classes
388 if((mapit->second).size()<minSize_opt[0]){
389 trainingMap.erase(mapit);
390 continue;
391 }
392 trainingPixels.push_back(mapit->second);
393 if(verbose_opt[0]>1)
394 std::cout << mapit->first << ": " << (mapit->second).size() << " samples" << std::endl;
395 ++mapit;
396 }
397 if(!ibag){
398 nclass=trainingPixels.size();
399 if(classname_opt.size())
400 assert(nclass==classname_opt.size());
401 nband=(training_opt.size())?trainingPixels[0][0].size()-2:trainingPixels[0][0].size();//X and Y
402 }
403 else{
404 assert(nclass==trainingPixels.size());
405 assert(nband==(training_opt.size())?trainingPixels[0][0].size()-2:trainingPixels[0][0].size());
406 }
407
408 //do not remove outliers here: could easily be obtained through ogr2ogr -where 'B2<110' output.shp input.shp
409 //balance training data
410 if(balance_opt[0]>0){
411 while(balance_opt.size()<nclass)
412 balance_opt.push_back(balance_opt.back());
413 if(random_opt[0])
414 srand(time(NULL));
415 totalSamples=0;
416 for(int iclass=0;iclass<nclass;++iclass){
417 if(trainingPixels[iclass].size()>balance_opt[iclass]){
418 while(trainingPixels[iclass].size()>balance_opt[iclass]){
419 int index=rand()%trainingPixels[iclass].size();
420 trainingPixels[iclass].erase(trainingPixels[iclass].begin()+index);
421 }
422 }
423 else{
424 int oldsize=trainingPixels[iclass].size();
425 for(int isample=trainingPixels[iclass].size();isample<balance_opt[iclass];++isample){
426 int index = rand()%oldsize;
427 trainingPixels[iclass].push_back(trainingPixels[iclass][index]);
428 }
429 }
430 totalSamples+=trainingPixels[iclass].size();
431 }
432 }
433
434 //set scale and offset
435 offset[ibag].resize(nband);
436 scale[ibag].resize(nband);
437 if(offset_opt.size()>1)
438 assert(offset_opt.size()==nband);
439 if(scale_opt.size()>1)
440 assert(scale_opt.size()==nband);
441 for(int iband=0;iband<nband;++iband){
442 if(verbose_opt[0]>=1)
443 cout << "scaling for band" << iband << endl;
444 offset[ibag][iband]=(offset_opt.size()==1)?offset_opt[0]:offset_opt[iband];
445 scale[ibag][iband]=(scale_opt.size()==1)?scale_opt[0]:scale_opt[iband];
446 //search for min and maximum
447 if(scale[ibag][iband]<=0){
448 float theMin=trainingPixels[0][0][iband+startBand];
449 float theMax=trainingPixels[0][0][iband+startBand];
450 for(int iclass=0;iclass<nclass;++iclass){
451 for(int isample=0;isample<trainingPixels[iclass].size();++isample){
452 if(theMin>trainingPixels[iclass][isample][iband+startBand])
453 theMin=trainingPixels[iclass][isample][iband+startBand];
454 if(theMax<trainingPixels[iclass][isample][iband+startBand])
455 theMax=trainingPixels[iclass][isample][iband+startBand];
456 }
457 }
458 offset[ibag][iband]=theMin+(theMax-theMin)/2.0;
459 scale[ibag][iband]=(theMax-theMin)/2.0;
460 if(verbose_opt[0]>=1){
461 std::cout << "Extreme image values for band " << iband << ": [" << theMin << "," << theMax << "]" << std::endl;
462 std::cout << "Using offset, scale: " << offset[ibag][iband] << ", " << scale[ibag][iband] << std::endl;
463 std::cout << "scaled values for band " << iband << ": [" << (theMin-offset[ibag][iband])/scale[ibag][iband] << "," << (theMax-offset[ibag][iband])/scale[ibag][iband] << "]" << std::endl;
464 }
465 }
466 }
467 }
468 else{//use same offset and scale
469 offset[ibag].resize(nband);
470 scale[ibag].resize(nband);
471 for(int iband=0;iband<nband;++iband){
472 offset[ibag][iband]=offset[0][iband];
473 scale[ibag][iband]=scale[0][iband];
474 }
475 }
476
477 if(!ibag){
478 if(priors_opt.size()==1){//default: equal priors for each class
479 priors.resize(nclass);
480 for(int iclass=0;iclass<nclass;++iclass)
481 priors[iclass]=1.0/nclass;
482 }
483 assert(priors_opt.size()==1||priors_opt.size()==nclass);
484
485 //set bagsize for each class if not done already via command line
486 while(bagSize_opt.size()<nclass)
487 bagSize_opt.push_back(bagSize_opt.back());
488
489 if(verbose_opt[0]>=1){
490 std::cout << "number of bands: " << nband << std::endl;
491 std::cout << "number of classes: " << nclass << std::endl;
492 std::cout << "priors:";
493 if(priorimg_opt.empty()){
494 for(int iclass=0;iclass<nclass;++iclass)
495 std::cout << " " << priors[iclass];
496 std::cout << std::endl;
497 }
498 }
499 map<string,Vector2d<float> >::iterator mapit=trainingMap.begin();
500 bool doSort=true;
501 while(mapit!=trainingMap.end()){
502 nameVector.push_back(mapit->first);
503 if(classValueMap.size()){
504 //check if name in training is covered by classname_opt (values can not be 0)
505 if(classValueMap[mapit->first]>0){
506 if(cm.getClassIndex(type2string<short>(classValueMap[mapit->first]))<0)
507 cm.pushBackClassName(type2string<short>(classValueMap[mapit->first]),doSort);
508 }
509 else{
510 std::cerr << "Error: names in classname option are not complete, please check names in training vector and make sure classvalue is > 0" << std::endl;
511 exit(1);
512 }
513 }
514 else
515 cm.pushBackClassName(mapit->first,doSort);
516 ++mapit;
517 }
518 if(classname_opt.empty()){
519 //std::cerr << "Warning: no class name and value pair provided for all " << nclass << " classes, using string2type<int> instead!" << std::endl;
520 for(int iclass=0;iclass<nclass;++iclass){
521 if(verbose_opt[0])
522 std::cout << iclass << " " << cm.getClass(iclass) << " -> " << string2type<short>(cm.getClass(iclass)) << std::endl;
523 classValueMap[cm.getClass(iclass)]=string2type<short>(cm.getClass(iclass));
524 }
525 }
526 if(priors_opt.size()==nameVector.size()){
527 std::cerr << "Warning: please check if priors are provided in correct order!!!" << std::endl;
528 for(int iclass=0;iclass<nameVector.size();++iclass)
529 std::cerr << nameVector[iclass] << " " << priors_opt[iclass] << std::endl;
530 }
531 }//if(!ibag)
532
533 //Calculate features of training set
534 vector< Vector2d<float> > trainingFeatures(nclass);
535 for(int iclass=0;iclass<nclass;++iclass){
536 int nctraining=0;
537 if(verbose_opt[0]>=1)
538 cout << "calculating features for class " << iclass << endl;
539 if(random_opt[0])
540 srand(time(NULL));
541 nctraining=(bagSize_opt[iclass]<100)? trainingPixels[iclass].size()/100.0*bagSize_opt[iclass] : trainingPixels[iclass].size();//bagSize_opt[iclass] given in % of training size
542 if(nctraining<=0)
543 nctraining=1;
544 assert(nctraining<=trainingPixels[iclass].size());
545 int index=0;
546 if(bagSize_opt[iclass]<100)
547 random_shuffle(trainingPixels[iclass].begin(),trainingPixels[iclass].end());
548
549 trainingFeatures[iclass].resize(nctraining);
550 for(int isample=0;isample<nctraining;++isample){
551 //scale pixel values according to scale and offset!!!
552 for(int iband=0;iband<nband;++iband){
553 float value=trainingPixels[iclass][isample][iband+startBand];
554 trainingFeatures[iclass][isample].push_back((value-offset[ibag][iband])/scale[ibag][iband]);
555 }
556 }
557 assert(trainingFeatures[iclass].size()==nctraining);
558 }
559
560 unsigned int nFeatures=trainingFeatures[0][0].size();
561 unsigned int ntraining=0;
562 for(int iclass=0;iclass<nclass;++iclass)
563 ntraining+=trainingFeatures[iclass].size();
564
565 const unsigned int num_layers = nneuron_opt.size()+2;
566 const float desired_error = 0.0003;
567 const unsigned int iterations_between_reports = (verbose_opt[0])? maxit_opt[0]+1:0;
568 if(verbose_opt[0]>=1){
569 cout << "number of features: " << nFeatures << endl;
570 cout << "creating artificial neural network with " << nneuron_opt.size() << " hidden layer, having " << endl;
571 for(int ilayer=0;ilayer<nneuron_opt.size();++ilayer)
572 cout << nneuron_opt[ilayer] << " ";
573 cout << "neurons" << endl;
574 cout << "connection_opt[0]: " << connection_opt[0] << std::endl;
575 cout << "num_layers: " << num_layers << std::endl;
576 cout << "nFeatures: " << nFeatures << std::endl;
577 cout << "nneuron_opt[0]: " << nneuron_opt[0] << std::endl;
578 cout << "number of classes (nclass): " << nclass << std::endl;
579 }
580 switch(num_layers){
581 case(3):{
582 // net[ibag].create_sparse(connection_opt[0],num_layers, nFeatures, nneuron_opt[0], nclass);//replace all create_sparse with create_sparse_array due to bug in FANN!
583 unsigned int layers[3];
584 layers[0]=nFeatures;
585 layers[1]=nneuron_opt[0];
586 layers[2]=nclass;
587 net[ibag].create_sparse_array(connection_opt[0],num_layers,layers);
588 break;
589 }
590 case(4):{
591 unsigned int layers[4];
592 layers[0]=nFeatures;
593 layers[1]=nneuron_opt[0];
594 layers[2]=nneuron_opt[1];
595 layers[3]=nclass;
596 // layers.push_back(nFeatures);
597 // for(int ihidden=0;ihidden<nneuron_opt.size();++ihidden)
598 // layers.push_back(nneuron_opt[ihidden]);
599 // layers.push_back(nclass);
600 net[ibag].create_sparse_array(connection_opt[0],num_layers,layers);
601 break;
602 }
603 default:
604 cerr << "Only 1 or 2 hidden layers are supported!" << endl;
605 exit(1);
606 break;
607 }
608 if(verbose_opt[0]>=1)
609 cout << "network created" << endl;
610
611 net[ibag].set_learning_rate(learning_opt[0]);
612
613 // net.set_activation_steepness_hidden(1.0);
614 // net.set_activation_steepness_output(1.0);
615
616 net[ibag].set_activation_function_hidden(FANN::SIGMOID_SYMMETRIC_STEPWISE);
617 net[ibag].set_activation_function_output(FANN::SIGMOID_SYMMETRIC_STEPWISE);
618
619 // Set additional properties such as the training algorithm
620 // net.set_training_algorithm(FANN::TRAIN_QUICKPROP);
621
622 // Output network type and parameters
623 if(verbose_opt[0]>=1){
624 cout << endl << "Network Type : ";
625 switch (net[ibag].get_network_type())
626 {
627 case FANN::LAYER:
628 cout << "LAYER" << endl;
629 break;
630 case FANN::SHORTCUT:
631 cout << "SHORTCUT" << endl;
632 break;
633 default:
634 cout << "UNKNOWN" << endl;
635 break;
636 }
637 net[ibag].print_parameters();
638 }
639
640 if(cv_opt[0]>1){
641 if(verbose_opt[0])
642 std::cout << "cross validation" << std::endl;
643 vector<unsigned short> referenceVector;
644 vector<unsigned short> outputVector;
645 float rmse=net[ibag].cross_validation(trainingFeatures,
646 ntraining,
647 cv_opt[0],
648 maxit_opt[0],
649 desired_error,
650 referenceVector,
651 outputVector,
652 verbose_opt[0]);
653 map<string,Vector2d<float> >::iterator mapit=trainingMap.begin();
654 for(int isample=0;isample<referenceVector.size();++isample){
655 string refClassName=nameVector[referenceVector[isample]];
656 string className=nameVector[outputVector[isample]];
657 if(classValueMap.size())
658 cm.incrementResult(type2string<short>(classValueMap[refClassName]),type2string<short>(classValueMap[className]),1.0/nbag);
659 else
660 cm.incrementResult(cm.getClass(referenceVector[isample]),cm.getClass(outputVector[isample]),1.0/nbag);
661 }
662 }
663
664 if(verbose_opt[0]>=1)
665 cout << endl << "Set training data" << endl;
666
667 if(verbose_opt[0]>=1)
668 cout << endl << "Training network" << endl;
669
670 if(verbose_opt[0]>=1){
671 cout << "Max Epochs " << setw(8) << maxit_opt[0] << ". "
672 << "Desired Error: " << left << desired_error << right << endl;
673 }
674 if(weights_opt.size()==net[ibag].get_total_connections()){//no new training needed (same training sample)
675 vector<fann_connection> convector;
676 net[ibag].get_connection_array(convector);
677 for(int i_connection=0;i_connection<net[ibag].get_total_connections();++i_connection)
678 convector[i_connection].weight=weights_opt[i_connection];
679 net[ibag].set_weight_array(convector);
680 }
681 else{
682 bool initWeights=true;
683 net[ibag].train_on_data(trainingFeatures,ntraining,initWeights, maxit_opt[0],
684 iterations_between_reports, desired_error);
685 }
686
687
688 if(verbose_opt[0]>=2){
689 net[ibag].print_connections();
690 vector<fann_connection> convector;
691 net[ibag].get_connection_array(convector);
692 for(int i_connection=0;i_connection<net[ibag].get_total_connections();++i_connection)
693 cout << "connection " << i_connection << ": " << convector[i_connection].weight << endl;
694
695 }
696 }//for ibag
697 if(cv_opt[0]>1){
698 assert(cm.nReference());
699 cm.setFormat(cmformat_opt[0]);
700 cm.reportSE95(false);
701 std::cout << cm << std::endl;
702 cout << "class #samples userAcc prodAcc" << endl;
703 double se95_ua=0;
704 double se95_pa=0;
705 double se95_oa=0;
706 double dua=0;
707 double dpa=0;
708 double doa=0;
709 for(int iclass=0;iclass<cm.nClasses();++iclass){
710 dua=cm.ua_pct(cm.getClass(iclass),&se95_ua);
711 dpa=cm.pa_pct(cm.getClass(iclass),&se95_pa);
712 cout << cm.getClass(iclass) << " " << cm.nReference(cm.getClass(iclass)) << " " << dua << " (" << se95_ua << ")" << " " << dpa << " (" << se95_pa << ")" << endl;
713 }
714 std::cout << "Kappa: " << cm.kappa() << std::endl;
715 doa=cm.oa_pct(&se95_oa);
716 std::cout << "Overall Accuracy: " << doa << " (" << se95_oa << ")" << std::endl;
717 }
718 //--------------------------------- end of training -----------------------------------
719 if(input_opt.empty())
720 exit(0);
721
722 const char* pszMessage;
723 void* pProgressArg=NULL;
724 GDALProgressFunc pfnProgress=GDALTermProgress;
725 float progress=0;
726 //-------------------------------- open image file ------------------------------------
727 bool inputIsRaster=false;
728 ImgReaderOgr imgReaderOgr;
729 try{
730 imgReaderOgr.open(input_opt[0]);
731 imgReaderOgr.close();
732 }
733 catch(string errorString){
734 inputIsRaster=true;
735 }
736 if(!inputIsRaster){//classify vector file
737 cm.clearResults();
738 //notice that fields have already been set by readDataImageOgr (taking into account appropriate bands)
739 for(int ivalidation=0;ivalidation<input_opt.size();++ivalidation){
740 if(output_opt.size())
741 assert(output_opt.size()==input_opt.size());
742 if(verbose_opt[0])
743 cout << "opening img reader " << input_opt[ivalidation] << endl;
744 imgReaderOgr.open(input_opt[ivalidation]);
745 ImgWriterOgr imgWriterOgr;
746
747 if(output_opt.size()){
748 if(verbose_opt[0])
749 std::cout << "opening img writer and copying fields from img reader" << output_opt[ivalidation] << std::endl;
750 imgWriterOgr.open(output_opt[ivalidation],imgReaderOgr);
751 }
752 if(verbose_opt[0])
753 cout << "number of layers in input ogr file: " << imgReaderOgr.getLayerCount() << endl;
754 for(int ilayer=0;ilayer<imgReaderOgr.getLayerCount();++ilayer){
755 if(verbose_opt[0])
756 cout << "processing input layer " << ilayer << endl;
757 if(output_opt.size()){
758 if(verbose_opt[0])
759 std::cout << "creating field class" << std::endl;
760 if(classValueMap.size())
761 imgWriterOgr.createField("class",OFTInteger,ilayer);
762 else
763 imgWriterOgr.createField("class",OFTString,ilayer);
764 }
765 unsigned int nFeatures=imgReaderOgr.getFeatureCount(ilayer);
766 unsigned int ifeature=0;
767 progress=0;
768 pfnProgress(progress,pszMessage,pProgressArg);
769 OGRFeature *poFeature;
770 while( (poFeature = imgReaderOgr.getLayer(ilayer)->GetNextFeature()) != NULL ){
771 if(verbose_opt[0]>1)
772 cout << "feature " << ifeature << endl;
773 if( poFeature == NULL ){
774 cout << "Warning: could not read feature " << ifeature << " in layer " << imgReaderOgr.getLayerName(ilayer) << endl;
775 continue;
776 }
777 OGRFeature *poDstFeature = NULL;
778 if(output_opt.size()){
779 poDstFeature=imgWriterOgr.createFeature(ilayer);
780 if( poDstFeature->SetFrom( poFeature, TRUE ) != OGRERR_NONE ){
781 CPLError( CE_Failure, CPLE_AppDefined,
782 "Unable to translate feature %d from layer %s.\n",
783 poFeature->GetFID(), imgWriterOgr.getLayerName(ilayer).c_str() );
784 OGRFeature::DestroyFeature( poFeature );
785 OGRFeature::DestroyFeature( poDstFeature );
786 }
787 }
788 vector<float> validationPixel;
789 vector<float> validationFeature;
790
791 imgReaderOgr.readData(validationPixel,OFTReal,fields,poFeature,ilayer);
792 assert(validationPixel.size()==nband);
793 vector<float> probOut(nclass);//posterior prob for each class
794 for(int iclass=0;iclass<nclass;++iclass)
795 probOut[iclass]=0;
796 for(int ibag=0;ibag<nbag;++ibag){
797 for(int iband=0;iband<nband;++iband){
798 validationFeature.push_back((validationPixel[iband]-offset[ibag][iband])/scale[ibag][iband]);
799 if(verbose_opt[0]==2)
800 std:: cout << " " << validationFeature.back();
801 }
802 if(verbose_opt[0]==2)
803 std::cout << std:: endl;
804 vector<float> result(nclass);
805 result=net[ibag].run(validationFeature);
806
807 if(verbose_opt[0]>1){
808 for(int iclass=0;iclass<result.size();++iclass)
809 std::cout << result[iclass] << " ";
810 std::cout << std::endl;
811 }
812 //calculate posterior prob of bag
813 for(int iclass=0;iclass<nclass;++iclass){
814 result[iclass]=(result[iclass]+1.0)/2.0;//bring back to scale [0,1]
815 switch(comb_opt[0]){
816 default:
817 case(0)://sum rule
818 probOut[iclass]+=result[iclass]*priors[iclass];//add probabilities for each bag
819 break;
820 case(1)://product rule
821 probOut[iclass]*=pow(static_cast<float>(priors[iclass]),static_cast<float>(1.0-nbag)/nbag)*result[iclass];//multiply probabilities for each bag
822 break;
823 case(2)://max rule
824 if(priors[iclass]*result[iclass]>probOut[iclass])
825 probOut[iclass]=priors[iclass]*result[iclass];
826 break;
827 }
828 }
829 }//for ibag
830 //search for max class prob
831 float maxBag=0;
832 float normBag=0;
833 string classOut="Unclassified";
834 for(int iclass=0;iclass<nclass;++iclass){
835 if(verbose_opt[0]>1)
836 std::cout << probOut[iclass] << " ";
837 if(probOut[iclass]>maxBag){
838 maxBag=probOut[iclass];
839 classOut=nameVector[iclass];
840 }
841 }
842 //look for class name
843 if(verbose_opt[0]>1){
844 if(classValueMap.size())
845 std::cout << "->" << classValueMap[classOut] << std::endl;
846 else
847 std::cout << "->" << classOut << std::endl;
848 }
849 if(output_opt.size()){
850 if(classValueMap.size())
851 poDstFeature->SetField("class",classValueMap[classOut]);
852 else
853 poDstFeature->SetField("class",classOut.c_str());
854 poDstFeature->SetFID( poFeature->GetFID() );
855 }
856 int labelIndex=poFeature->GetFieldIndex(label_opt[0].c_str());
857 if(labelIndex>=0){
858 string classRef=poFeature->GetFieldAsString(labelIndex);
859 if(classRef!="0"){
860 if(classValueMap.size())
861 cm.incrementResult(type2string<short>(classValueMap[classRef]),type2string<short>(classValueMap[classOut]),1);
862 else
863 cm.incrementResult(classRef,classOut,1);
864 }
865 }
866 CPLErrorReset();
867 if(output_opt.size()){
868 if(imgWriterOgr.createFeature(poDstFeature,ilayer) != OGRERR_NONE){
869 CPLError( CE_Failure, CPLE_AppDefined,
870 "Unable to translate feature %d from layer %s.\n",
871 poFeature->GetFID(), imgWriterOgr.getLayerName(ilayer).c_str() );
872 OGRFeature::DestroyFeature( poDstFeature );
873 OGRFeature::DestroyFeature( poDstFeature );
874 }
875 }
876 ++ifeature;
877 if(!verbose_opt[0]){
878 progress=static_cast<float>(ifeature+1.0)/nFeatures;
879 pfnProgress(progress,pszMessage,pProgressArg);
880 }
881 OGRFeature::DestroyFeature( poFeature );
882 OGRFeature::DestroyFeature( poDstFeature );
883 }//get next feature
884 }//next layer
885 imgReaderOgr.close();
886 if(output_opt.size())
887 imgWriterOgr.close();
888 }
889 if(cm.nReference()){
890 std::cout << cm << std::endl;
891 // cout << "class #samples userAcc prodAcc" << endl;
892 // double se95_ua=0;
893 // double se95_pa=0;
894 // double se95_oa=0;
895 // double dua=0;
896 // double dpa=0;
897 // double doa=0;
898 // for(short iclass=0;iclass<cm.nClasses();++iclass){
899 // dua=cm.ua_pct(cm.getClass(iclass),&se95_ua);
900 // dpa=cm.pa_pct(cm.getClass(iclass),&se95_pa);
901 // cout << cm.getClass(iclass) << " " << cm.nReference(cm.getClass(iclass)) << " " << dua << " (" << se95_ua << ")" << " " << dpa << " (" << se95_pa << ")" << endl;
902 // }
903 // std::cout << "Kappa: " << cm.kappa() << std::endl;
904 // doa=cm.oa_pct(&se95_oa);
905 // std::cout << "Overall Accuracy: " << doa << " (" << se95_oa << ")" << std::endl;
906 }
907 }
908 try{
909 if(active_opt.size())
910 activeWriter.close();
911 }
912 catch(string errorString){
913 std::cerr << "Error: errorString" << std::endl;
914 }
915 return 0;
916}