pktools  2.6.7
Processing Kernel for geospatial data
pkdiff.cc
1 /**********************************************************************
2 pkdiff.cc: program to compare two raster image files
3 Copyright (C) 2008-2014 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 <assert.h>
21 #include "imageclasses/ImgReaderGdal.h"
22 #include "imageclasses/ImgWriterGdal.h"
23 #include "imageclasses/ImgReaderOgr.h"
24 #include "imageclasses/ImgWriterOgr.h"
25 #include "base/Optionpk.h"
26 #include "algorithms/ConfusionMatrix.h"
27 
28 /******************************************************************************/
93 using namespace std;
94 
95 int main(int argc, char *argv[])
96 {
97  Optionpk<string> input_opt("i", "input", "Input raster dataset.");
98  Optionpk<string> reference_opt("ref", "reference", "Reference (raster or vector) dataset");
99  Optionpk<string> layer_opt("ln", "ln", "Layer name(s) in sample. Leave empty to select all (for vector reference datasets only)");
100  Optionpk<string> mask_opt("m", "mask", "Use the first band of the specified file as a validity mask. Nodata values can be set with the option msknodata.");
101  Optionpk<double> msknodata_opt("msknodata", "msknodata", "Mask value(s) where image is invalid. Use negative value for valid data (example: use -t -1: if only -1 is valid value)", 0);
102  Optionpk<double> nodata_opt("nodata", "nodata", "No data value(s) in input or reference dataset are ignored");
103  Optionpk<short> band_opt("b", "band", "Input (reference) raster band. Optionally, you can define different bands for input and reference bands respectively: -b 1 -b 0.", 0);
104  Optionpk<bool> rmse_opt("rmse", "rmse", "Report root mean squared error", false);
105  Optionpk<bool> regression_opt("reg", "reg", "Report linear regression (Input = c0+c1*Reference)", false);
106  Optionpk<bool> confusion_opt("cm", "confusion", "Create confusion matrix (to std out)", false);
107  Optionpk<string> cmformat_opt("cmf","cmf","Format for confusion matrix (ascii or latex)","ascii");
108  Optionpk<string> cmoutput_opt("cmo","cmo","Output file for confusion matrix");
109  Optionpk<bool> se95_opt("se95","se95","Report standard error for 95 confidence interval",false);
110  Optionpk<string> labelref_opt("lr", "lref", "Attribute name of the reference label (for vector reference datasets only)", "label");
111  Optionpk<string> classname_opt("c", "class", "List of class names.");
112  Optionpk<short> classvalue_opt("r", "reclass", "List of class values (use same order as in classname option).");
113  Optionpk<string> output_opt("o", "output", "Output dataset (optional)");
114  Optionpk<string> ogrformat_opt("f", "f", "OGR format for output vector (for vector reference datasets only)","SQLite");
115  Optionpk<string> oformat_opt("of", "oformat", "Output image format (see also gdal_translate).","GTiff");
116  Optionpk<string> labelclass_opt("lc", "lclass", "Attribute name of the classified label (for vector reference datasets only)", "class");
117  Optionpk<short> boundary_opt("bnd", "boundary", "Boundary for selecting the sample (for vector reference datasets only)", 1,1);
118  Optionpk<bool> homogeneous_opt("hom", "homogeneous", "Only take regions with homogeneous boundary into account (for reference datasets only)", false,1);
119  Optionpk<bool> disc_opt("circ", "circular", "Use circular boundary (for vector reference datasets only)", false,1);
120  Optionpk<string> colorTable_opt("ct", "ct", "Color table in ASCII format having 5 columns: id R G B ALFA (0: transparent, 255: solid).");
121  Optionpk<string> option_opt("co", "co", "Creation option for output file. Multiple options can be specified.");
122  Optionpk<short> valueE_opt("\0", "correct", "Value for correct pixels", 0,2);
123  Optionpk<short> valueO_opt("\0", "omission", "Value for omission errors: input label > reference label", 1,2);
124  Optionpk<short> valueC_opt("\0", "commission", "Value for commission errors: input label < reference label", 2,1);
125  Optionpk<short> verbose_opt("v", "verbose", "Verbose level", 0,2);
126 
127  output_opt.setHide(1);
128  ogrformat_opt.setHide(1);
129  oformat_opt.setHide(1);
130  labelclass_opt.setHide(1);
131  boundary_opt.setHide(1);
132  homogeneous_opt.setHide(1);
133  disc_opt.setHide(1);
134  colorTable_opt.setHide(1);
135  option_opt.setHide(1);
136 
137  bool doProcess;//stop process when program was invoked with help option (-h --help)
138  try{
139  doProcess=input_opt.retrieveOption(argc,argv);
140  reference_opt.retrieveOption(argc,argv);
141  layer_opt.retrieveOption(argc,argv);
142  band_opt.retrieveOption(argc,argv);
143  rmse_opt.retrieveOption(argc,argv);
144  regression_opt.retrieveOption(argc,argv);
145  confusion_opt.retrieveOption(argc,argv);
146  labelref_opt.retrieveOption(argc,argv);
147  classname_opt.retrieveOption(argc,argv);
148  classvalue_opt.retrieveOption(argc,argv);
149  nodata_opt.retrieveOption(argc,argv);
150  mask_opt.retrieveOption(argc,argv);
151  msknodata_opt.retrieveOption(argc,argv);
152  output_opt.retrieveOption(argc,argv);
153  ogrformat_opt.retrieveOption(argc,argv);
154  labelclass_opt.retrieveOption(argc,argv);
155  cmformat_opt.retrieveOption(argc,argv);
156  cmoutput_opt.retrieveOption(argc,argv);
157  se95_opt.retrieveOption(argc,argv);
158  boundary_opt.retrieveOption(argc,argv);
159  homogeneous_opt.retrieveOption(argc,argv);
160  disc_opt.retrieveOption(argc,argv);
161  colorTable_opt.retrieveOption(argc,argv);
162  option_opt.retrieveOption(argc,argv);
163  // class_opt.retrieveOption(argc,argv);
164  valueE_opt.retrieveOption(argc,argv);
165  valueO_opt.retrieveOption(argc,argv);
166  valueC_opt.retrieveOption(argc,argv);
167  verbose_opt.retrieveOption(argc,argv);
168  }
169  catch(string predefinedString){
170  std::cout << predefinedString << std::endl;
171  exit(0);
172  }
173  if(!doProcess){
174  cout << endl;
175  cout << "Usage: pkdiff -i input -ref reference" << endl;
176  cout << endl;
177  std::cout << "short option -h shows basic options only, use long option --help to show all options" << std::endl;
178  exit(0);//help was invoked, stop processing
179  }
180 
181  ImgReaderGdal inputReader;
182  ImgReaderGdal maskReader;
183 
184  if(verbose_opt[0]){
185  cout << "flag(s) set to";
186  for(int iflag=0;iflag<nodata_opt.size();++iflag)
187  cout << " " << nodata_opt[iflag];
188  cout << endl;
189  }
190 
191  if(input_opt.empty()){
192  std::cerr << "No input file provided (use option -i). Use --help for help information" << std::endl;
193  exit(0);
194  }
195  if(reference_opt.empty()){
196  std::cerr << "No reference file provided (use option -ref). Use --help for help information" << std::endl;
197  exit(0);
198  }
199 
200  //band_opt[0] is for input
201  //band_opt[1] is for reference
202  if(band_opt.size()<2)
203  band_opt.push_back(band_opt[0]);
204 
205  if(mask_opt.size())
206  while(mask_opt.size()<input_opt.size())
207  mask_opt.push_back(mask_opt[0]);
208  vector<short> inputRange;
209  vector<short> referenceRange;
211  int nclass=0;
212  map<string,short> classValueMap;
213  vector<std::string> nameVector(255);//the inverse of the classValueMap
214  vector<string> classNames;
215 
216  unsigned int ntotalValidation=0;
217  unsigned int nflagged=0;
218  Vector2d<int> resultClass;
219  vector<float> user;
220  vector<float> producer;
221  vector<unsigned int> nvalidation;
222 
223  if(confusion_opt[0]){
224  // if(class_opt.size()>1)
225  // inputRange=class_opt;
226  // if(classvalue_opt.size()>1)
227  // inputRange=classvalue_opt;
228  // else{
229  try{
230  if(verbose_opt[0])
231  cout << "opening input image file " << input_opt[0] << endl;
232  inputReader.open(input_opt[0]);//,imagicX_opt[0],imagicY_opt[0]);
233  }
234  catch(string error){
235  cerr << error << endl;
236  exit(1);
237  }
238  inputReader.getRange(inputRange,band_opt[0]);
239  inputReader.close();
240  // }
241 
242  for(int iflag=0;iflag<nodata_opt.size();++iflag){
243  vector<short>::iterator fit;
244  fit=find(inputRange.begin(),inputRange.end(),static_cast<short>(nodata_opt[iflag]));
245  if(fit!=inputRange.end())
246  inputRange.erase(fit);
247  }
248  nclass=inputRange.size();
249  if(verbose_opt[0]){
250  cout << "nclass (inputRange.size()): " << nclass << endl;
251  cout << "input range: " << endl;
252  }
253  if(classname_opt.size()){
254  assert(classname_opt.size()==classvalue_opt.size());
255  for(int iclass=0;iclass<classname_opt.size();++iclass){
256  classValueMap[classname_opt[iclass]]=classvalue_opt[iclass];
257  assert(classvalue_opt[iclass]<nameVector.size());
258  nameVector[classvalue_opt[iclass]]=classname_opt[iclass];
259  }
260  }
261  // nclass=classValueMap.size();
262  for(int rc=0;rc<inputRange.size();++rc){
263  classNames.push_back(type2string(inputRange[rc]));
264  if(verbose_opt[0])
265  cout << inputRange[rc] << endl;
266  }
267  cm.setClassNames(classNames);
268  if(verbose_opt[0]){
269  cout << "class names: " << endl;
270  for(int iclass=0;iclass<cm.nClasses();++iclass)
271  cout << iclass << " " << cm.getClass(iclass) << endl;
272  }
273  resultClass.resize(nclass,nclass);
274  user.resize(nclass);
275  producer.resize(nclass);
276  nvalidation.resize(nclass);
277  //initialize
278  for(int rc=0;rc<nclass;++rc){
279  for(int ic=0;ic<nclass;++ic)
280  resultClass[rc][ic]=0;
281  nvalidation[rc]=0;
282  }
283  }
284 
285  bool isDifferent=false;
286  bool refIsRaster=false;
287 
288  ImgReaderOgr referenceReaderOgr;
289  try{
290  referenceReaderOgr.open(reference_opt[0]);
291  double ulx,uly,lrx,lry;
292  if(!referenceReaderOgr.getExtent(ulx,uly,lrx,lry))
293  throw(std::string("Input is raster"));
294  referenceReaderOgr.close();
295  if(verbose_opt[0])
296  std::cout << "reference is vector dataset" << std::endl;
297  }
298  catch(string errorString){
299  //todo: sampleIsRaster will not work from GDAL 2.0!!?? (unification of driver for raster and vector datasets)
300  refIsRaster=true;
301  if(verbose_opt[0])
302  std::cout << "reference is raster dataset" << std::endl;
303  }
304  const char* pszMessage;
305  void* pProgressArg=NULL;
306  GDALProgressFunc pfnProgress=GDALTermProgress;
307  float progress=0;
308  // if(reference_opt[0].find(".shp")!=string::npos){
309  if(!refIsRaster){
310  for(int iinput=0;iinput<input_opt.size();++iinput){
311  if(verbose_opt[0])
312  cout << "Processing input " << input_opt[iinput] << endl;
313  if(output_opt.size())
314  assert(reference_opt.size()==output_opt.size());
315  for(int iref=0;iref<reference_opt.size();++iref){
316  cout << "reference " << reference_opt[iref] << endl;
317  // assert(reference_opt[iref].find(".shp")!=string::npos);
318  try{
319  inputReader.open(input_opt[iinput]);//,imagicX_opt[0],imagicY_opt[0]);
320  if(mask_opt.size()){
321  maskReader.open(mask_opt[iinput]);
322  assert(inputReader.nrOfCol()==maskReader.nrOfCol());
323  assert(inputReader.nrOfRow()==maskReader.nrOfRow());
324  }
325  referenceReaderOgr.open(reference_opt[iref]);
326  }
327  catch(string error){
328  cerr << error << endl;
329  exit(1);
330  }
331  if(confusion_opt[0])
332  referenceRange=inputRange;
333 
334  ImgWriterOgr ogrWriter;
335  if(output_opt.size()){
336  try{
337  ogrWriter.open(output_opt[iref],ogrformat_opt[0]);
338  }
339  catch(string error){
340  cerr << error << endl;
341  exit(1);
342  }
343  }
344  int nlayer=referenceReaderOgr.getDataSource()->GetLayerCount();
345  for(int ilayer=0;ilayer<nlayer;++ilayer){
346  progress=0;
347  OGRLayer *readLayer=referenceReaderOgr.getLayer(ilayer);
348  // readLayer = referenceReaderOgr.getDataSource()->GetLayer(ilayer);
349  string currentLayername=readLayer->GetName();
350  if(layer_opt.size())
351  if(find(layer_opt.begin(),layer_opt.end(),currentLayername)==layer_opt.end())
352  continue;
353  if(!verbose_opt[0])
354  pfnProgress(progress,pszMessage,pProgressArg);
355  else
356  cout << "processing layer " << readLayer->GetName() << endl;
357 
358  readLayer->ResetReading();
359  OGRLayer *writeLayer;
360  if(output_opt.size()){
361  if(verbose_opt[0])
362  cout << "creating output vector file " << output_opt[0] << endl;
363  // assert(output_opt[0].find(".shp")!=string::npos);
364  char **papszOptions=NULL;
365  if(verbose_opt[0])
366  cout << "creating layer: " << readLayer->GetName() << endl;
367  // if(ogrWriter.createLayer(layername, referenceReaderOgr.getProjection(ilayer), referenceReaderOgr.getGeometryType(ilayer), papszOptions)==NULL)
368  writeLayer=ogrWriter.createLayer(readLayer->GetName(), referenceReaderOgr.getProjection(ilayer), wkbPoint, papszOptions);
369  assert(writeLayer);
370  if(verbose_opt[0]){
371  cout << "created layer" << endl;
372  cout << "copy fields from " << reference_opt[iref] << endl;
373  }
374  ogrWriter.copyFields(referenceReaderOgr,ilayer,ilayer);
375  //create extra field for classified label
376  short theDim=boundary_opt[0];
377  for(int windowJ=-theDim/2;windowJ<(theDim+1)/2;++windowJ){
378  for(int windowI=-theDim/2;windowI<(theDim+1)/2;++windowI){
379  if(disc_opt[0]&&(windowI*windowI+windowJ*windowJ>(theDim/2)*(theDim/2)))
380  continue;
381  ostringstream fs;
382  if(theDim>1)
383  fs << labelclass_opt[0] << "_" << windowJ << "_" << windowI;
384  else
385  fs << labelclass_opt[0];
386  if(verbose_opt[0])
387  cout << "creating field " << fs.str() << endl;
388  ogrWriter.createField(fs.str(),OFTInteger,ilayer);
389  }
390  }
391  }
392  OGRFeature *readFeature;
393  OGRFeature *writeFeature;
394  int isample=0;
395  unsigned int nfeatureInLayer=readLayer->GetFeatureCount();
396  unsigned int ifeature=0;
397  while( (readFeature = readLayer->GetNextFeature()) != NULL ){
398  if(verbose_opt[0])
399  cout << "sample " << ++isample << endl;
400  //get x and y from readFeature
401  double x,y;
402  OGRGeometry *poGeometry;
403  OGRPoint centroidPoint;
404  OGRPoint *poPoint;
405  poGeometry = readFeature->GetGeometryRef();
406  // assert( poGeometry != NULL && wkbFlatten(poGeometry->getGeometryType()) == wkbPoint );
407  if(poGeometry==NULL)
408  continue;
409  else if(wkbFlatten(poGeometry->getGeometryType()) == wkbMultiPolygon){
410  OGRMultiPolygon readPolygon = *((OGRMultiPolygon *) poGeometry);
411  readPolygon = *((OGRMultiPolygon *) poGeometry);
412  readPolygon.Centroid(&centroidPoint);
413  poPoint=&centroidPoint;
414  }
415  else if(wkbFlatten(poGeometry->getGeometryType()) == wkbPolygon){
416  OGRPolygon readPolygon=*((OGRPolygon *) poGeometry);
417  readPolygon.Centroid(&centroidPoint);
418  poPoint=&centroidPoint;
419  }
420  else if(wkbFlatten(poGeometry->getGeometryType()) == wkbPoint )
421  poPoint = (OGRPoint *) poGeometry;
422  else{
423  std::cerr << "Warning: skipping feature (not of type point or polygon)" << std::endl;
424  continue;
425  }
426  x=poPoint->getX();
427  y=poPoint->getY();
428  double inputValue;
429  vector<double> inputValues;
430  bool isHomogeneous=true;
431  short maskValue;
432  short outputValue;
433  //read referenceValue from feature
434  unsigned short referenceValue;
435  string referenceClassName;
436  if(classValueMap.size()){
437  referenceClassName=readFeature->GetFieldAsString(readFeature->GetFieldIndex(labelref_opt[0].c_str()));
438  referenceValue=classValueMap[referenceClassName];
439  }
440  else
441  referenceValue=readFeature->GetFieldAsInteger(readFeature->GetFieldIndex(labelref_opt[0].c_str()));
442  if(verbose_opt[0])
443  cout << "reference value: " << referenceValue << endl;
444 
445  bool pixelFlagged=false;
446  bool maskFlagged=false;
447  for(int iflag=0;iflag<nodata_opt.size();++iflag){
448  if(referenceValue==nodata_opt[iflag])
449  pixelFlagged=true;
450  }
451  if(pixelFlagged)
452  continue;
453  double i_centre,j_centre;
454  //input reader is georeferenced!
455  inputReader.geo2image(x,y,i_centre,j_centre);
456  // else{
457  // i_centre=x;
458  // j_centre=y;
459  // }
460  //nearest neighbour
461  j_centre=static_cast<int>(j_centre);
462  i_centre=static_cast<int>(i_centre);
463  //check if j_centre is out of bounds
464  if(static_cast<int>(j_centre)<0||static_cast<int>(j_centre)>=inputReader.nrOfRow())
465  continue;
466  //check if i_centre is out of bounds
467  if(static_cast<int>(i_centre)<0||static_cast<int>(i_centre)>=inputReader.nrOfCol())
468  continue;
469 
470  if(output_opt.size()){
471  writeFeature = OGRFeature::CreateFeature(writeLayer->GetLayerDefn());
472  assert(readFeature);
473  int nfield=readFeature->GetFieldCount();
474  writeFeature->SetGeometry(poPoint);
475  if(verbose_opt[0])
476  cout << "copying fields from " << reference_opt[0] << endl;
477  assert(readFeature);
478  assert(writeFeature);
479  vector<int> panMap(nfield);
480  vector<int>::iterator panit=panMap.begin();
481  for(int ifield=0;ifield<nfield;++ifield)
482  panMap[ifield]=ifield;
483  writeFeature->SetFieldsFrom(readFeature,&(panMap[0]));
484  // if(writeFeature->SetFrom(readFeature)!= OGRERR_NONE)
485  // cerr << "writing feature failed" << endl;
486  // if(verbose_opt[0])
487  // cout << "feature written" << endl;
488  }
489  bool windowAllFlagged=true;
490  bool windowHasFlag=false;
491  short theDim=boundary_opt[0];
492  for(int windowJ=-theDim/2;windowJ<(theDim+1)/2;++windowJ){
493  for(int windowI=-theDim/2;windowI<(theDim+1)/2;++windowI){
494  if(disc_opt[0]&&(windowI*windowI+windowJ*windowJ>(theDim/2)*(theDim/2)))
495  continue;
496  int j=j_centre+windowJ;
497  //check if j is out of bounds
498  if(static_cast<int>(j)<0||static_cast<int>(j)>=inputReader.nrOfRow())
499  continue;
500  int i=i_centre+windowI;
501  //check if i is out of bounds
502  if(static_cast<int>(i)<0||static_cast<int>(i)>=inputReader.nrOfCol())
503  continue;
504  if(verbose_opt[0])
505  cout << setprecision(12) << "reading image value at x,y " << x << "," << y << " (" << i << "," << j << "), ";
506  inputReader.readData(inputValue,i,j,band_opt[0]);
507  inputValues.push_back(inputValue);
508  if(inputValues.back()!=*(inputValues.begin()))
509  isHomogeneous=false;
510  if(verbose_opt[0])
511  cout << "input value: " << inputValue << endl;
512  pixelFlagged=false;
513  for(int iflag=0;iflag<nodata_opt.size();++iflag){
514  if(inputValue==nodata_opt[iflag]){
515  pixelFlagged=true;
516  break;
517  }
518  }
519  maskFlagged=false;//(msknodata_opt[ivalue]>=0)?false:true;
520  if(mask_opt.size()){
521  maskReader.readData(maskValue,i,j,0);
522  for(int ivalue=0;ivalue<msknodata_opt.size();++ivalue){
523  if(msknodata_opt[ivalue]>=0){//values set in msknodata_opt are invalid
524  if(maskValue==msknodata_opt[ivalue]){
525  maskFlagged=true;
526  break;
527  }
528  }
529  else{//only values set in msknodata_opt are valid
530  if(maskValue!=-msknodata_opt[ivalue])
531  maskFlagged=true;
532  else{
533  maskFlagged=false;
534  break;
535  }
536  }
537  }
538  }
539  pixelFlagged=pixelFlagged||maskFlagged;
540  if(pixelFlagged)
541  windowHasFlag=true;
542  else
543  windowAllFlagged=false;//at least one good pixel in neighborhood
544  }
545  }
546  //at this point we know the values for the entire window
547 
548  if(homogeneous_opt[0]){//only centre pixel
549  int j=j_centre;
550  int i=i_centre;
551  //flag if not all pixels are homogeneous or if at least one pixel flagged
552 
553  if(!windowHasFlag&&isHomogeneous){
554  if(output_opt.size())
555  writeFeature->SetField(labelclass_opt[0].c_str(),static_cast<int>(inputValue));
556  if(confusion_opt[0]){
557  ++ntotalValidation;
558  if(classValueMap.size()){
559  assert(inputValue<nameVector.size());
560  string className=nameVector[static_cast<unsigned short>(inputValue)];
561  cm.incrementResult(type2string<short>(classValueMap[referenceClassName]),type2string<short>(classValueMap[className]),1);
562  }
563  else{
564  int rc=distance(referenceRange.begin(),find(referenceRange.begin(),referenceRange.end(),static_cast<unsigned short>(referenceValue)));
565  int ic=distance(inputRange.begin(),find(inputRange.begin(),inputRange.end(),static_cast<unsigned short>(inputValue)));
566  assert(rc<nclass);
567  assert(ic<nclass);
568  ++nvalidation[rc];
569  ++resultClass[rc][ic];
570  if(verbose_opt[0]>1)
571  cout << "increment: " << rc << " " << referenceRange[rc] << " " << ic << " " << inputRange[ic] << endl;
572  cm.incrementResult(cm.getClass(rc),cm.getClass(ic),1);
573  }
574  }
575  if(inputValue==referenceValue){//correct
576  outputValue=valueE_opt[0];
577  if(nodata_opt.size()){
578  if(valueE_opt[0]==nodata_opt[0])
579  outputValue=inputValue;
580  }
581  }
582  else if(inputValue>referenceValue)//1=forest,2=non-forest
583  outputValue=valueO_opt[0];//omission error
584  else
585  outputValue=valueC_opt[0];//commission error
586  }
587  }
588  else{
589  for(int windowJ=-theDim/2;windowJ<(theDim+1)/2;++windowJ){
590  for(int windowI=-theDim/2;windowI<(theDim+1)/2;++windowI){
591  if(disc_opt[0]&&(windowI*windowI+windowJ*windowJ>(theDim/2)*(theDim/2)))
592  continue;
593  int j=j_centre+windowJ;
594  //check if j is out of bounds
595  if(static_cast<int>(j)<0||static_cast<int>(j)>=inputReader.nrOfRow())
596  continue;
597  int i=i_centre+windowI;
598  //check if i is out of bounds
599  if(static_cast<int>(i)<0||static_cast<int>(i)>=inputReader.nrOfCol())
600  continue;
601  if(!windowAllFlagged){
602  ostringstream fs;
603  if(theDim>1)
604  fs << labelclass_opt[0] << "_" << windowJ << "_" << windowI;
605  else
606  fs << labelclass_opt[0];
607  if(output_opt.size())
608  writeFeature->SetField(fs.str().c_str(),inputValue);
609  if(!windowJ&&!windowI){//centre pixel
610  if(confusion_opt[0]){
611  ++ntotalValidation;
612  if(classValueMap.size()){
613  assert(inputValue<nameVector.size());
614  string className=nameVector[static_cast<unsigned short>(inputValue)];
615  cm.incrementResult(type2string<short>(classValueMap[referenceClassName]),type2string<short>(classValueMap[className]),1);
616  }
617  else{
618  int rc=distance(referenceRange.begin(),find(referenceRange.begin(),referenceRange.end(),static_cast<unsigned short>(referenceValue)));
619  int ic=distance(inputRange.begin(),find(inputRange.begin(),inputRange.end(),static_cast<unsigned short>(inputValue)));
620  if(rc>=nclass)
621  continue;
622  if(ic>=nclass)
623  continue;
624  // assert(rc<nclass);
625  // assert(ic<nclass);
626  ++nvalidation[rc];
627  ++resultClass[rc][ic];
628  if(verbose_opt[0]>1)
629  cout << "increment: " << rc << " " << referenceRange[rc] << " " << ic << " " << inputRange[ic] << endl;
630  cm.incrementResult(cm.getClass(rc),cm.getClass(ic),1);
631  }
632  }
633  if(inputValue==referenceValue){//correct
634  outputValue=valueE_opt[0];
635  if(nodata_opt.size()){
636  if(valueE_opt[0]==nodata_opt[0])
637  outputValue=inputValue;
638  }
639  }
640  else if(inputValue>referenceValue)//1=forest,2=non-forest
641  outputValue=valueO_opt[0];//omission error
642  else
643  outputValue=valueC_opt[0];//commission error
644  }
645  }
646  }
647  }
648  }
649  if(output_opt.size()){
650  if(!windowAllFlagged){
651  if(verbose_opt[0])
652  cout << "creating feature" << endl;
653  if(writeLayer->CreateFeature( writeFeature ) != OGRERR_NONE ){
654  string errorString="Failed to create feature in OGR vector file";
655  throw(errorString);
656  }
657  }
658  OGRFeature::DestroyFeature( writeFeature );
659  }
660  ++ifeature;
661  progress=static_cast<float>(ifeature+1)/nfeatureInLayer;
662  pfnProgress(progress,pszMessage,pProgressArg);
663  }//next feature
664  }//next layer
665  if(output_opt.size())
666  ogrWriter.close();
667  referenceReaderOgr.close();
668  inputReader.close();
669  if(mask_opt.size())
670  maskReader.close();
671  }//next reference
672  }//next input
673  pfnProgress(1.0,pszMessage,pProgressArg);
674  }//reference is OGR vector
675  else{//reference is GDAL raster
676  ImgWriterGdal gdalWriter;
677  try{
678  inputReader.open(input_opt[0]);
679  if(mask_opt.size())
680  maskReader.open(mask_opt[0]);
681  if(output_opt.size()){
682  if(verbose_opt[0])
683  cout << "opening output image " << output_opt[0] << endl;
684  if(option_opt.findSubstring("INTERLEAVE=")==option_opt.end()){
685  string theInterleave="INTERLEAVE=";
686  theInterleave+=inputReader.getInterleave();
687  option_opt.push_back(theInterleave);
688  }
689  gdalWriter.open(output_opt[0],inputReader.nrOfCol(),inputReader.nrOfRow(),1,inputReader.getDataType(),oformat_opt[0],option_opt);
690  if(nodata_opt.size())
691  gdalWriter.GDALSetNoDataValue(nodata_opt[0]);
692  gdalWriter.copyGeoTransform(inputReader);
693  if(colorTable_opt.size())
694  gdalWriter.setColorTable(colorTable_opt[0]);
695  else if(inputReader.getColorTable()!=NULL){
696  if(verbose_opt[0])
697  cout << "set colortable from input image" << endl;
698  gdalWriter.setColorTable(inputReader.getColorTable());
699  }
700  }
701  else if(verbose_opt[0])
702  cout << "no output image defined" << endl;
703 
704  }
705  catch(string error){
706  cout << error << endl;
707  exit(2);
708  }
709  //todo: support different data types!
710  vector<double> lineInput(inputReader.nrOfCol());
711  vector<double> lineMask(maskReader.nrOfCol());
712  vector<double> lineOutput;
713  vector<double> bufferInput;//for regression
714  vector<double> bufferReference;//for regression
715  if(output_opt.size())
716  lineOutput.resize(inputReader.nrOfCol());
717 
718  int irow=0;
719  int icol=0;
720  double oldreferencerow=-1;
721  double oldmaskrow=-1;
722  ImgReaderGdal referenceReaderGdal;
723  try{
724  referenceReaderGdal.open(reference_opt[0]);//,rmagicX_opt[0],rmagicY_opt[0]);
725  }
726  catch(string error){
727  cerr << error << endl;
728  exit(1);
729  }
730  if(inputReader.isGeoRef()){
731  assert(referenceReaderGdal.isGeoRef());
732  if(inputReader.getProjection()!=referenceReaderGdal.getProjection())
733  cerr << "Warning: projection of input image and reference image are different" << endl;
734  }
735  vector<double> lineReference(referenceReaderGdal.nrOfCol());
736  if(confusion_opt[0]){
737  referenceReaderGdal.getRange(referenceRange,band_opt[1]);
738  for(int iflag=0;iflag<nodata_opt.size();++iflag){
739  vector<short>::iterator fit;
740  fit=find(referenceRange.begin(),referenceRange.end(),static_cast<unsigned short>(nodata_opt[iflag]));
741  if(fit!=referenceRange.end())
742  referenceRange.erase(fit);
743  }
744  if(verbose_opt[0]){
745  cout << "reference range: " << endl;
746  for(int rc=0;rc<referenceRange.size();++rc)
747  cout << referenceRange[rc] << endl;
748  }
749  if(referenceRange.size()!=inputRange.size()){
750  if(confusion_opt[0]||output_opt.size()){
751  cout << "reference range is not equal to input range!" << endl;
752  cout << "Kappa: " << 0 << endl;
753  cout << "total weighted: " << 0 << endl;
754  }
755  else
756  cout << "reference range is not equal to input range!" << endl;
757  cout << input_opt[0] << " and " << reference_opt[0] << " are different" << endl;
758  exit(1);
759  }
760  }
761  double rmse=0;
762  // for(irow=0;irow<inputReader.nrOfRow()&&!isDifferent;++irow){
763  for(irow=0;irow<inputReader.nrOfRow();++irow){
764  //read line in lineInput, lineReference and lineMask
765  inputReader.readData(lineInput,irow,band_opt[0]);
766  double x,y;//geo coordinates
767  double ireference,jreference;//image coordinates in reference image
768  double imask,jmask;//image coordinates in mask image
769  for(icol=0;icol<inputReader.nrOfCol();++icol){
770  //find col in reference
771  inputReader.image2geo(icol,irow,x,y);
772  referenceReaderGdal.geo2image(x,y,ireference,jreference);
773  if(ireference<0||ireference>=referenceReaderGdal.nrOfCol()){
774  if(rmse_opt[0]||regression_opt[0])
775  continue;
776  else{
777  cerr << ireference << " out of reference range!" << endl;
778  cerr << x << " " << y << " " << icol << " " << irow << endl;
779  cerr << x << " " << y << " " << ireference << " " << jreference << endl;
780  exit(1);
781  }
782  }
783  if(jreference!=oldreferencerow){
784  if(jreference<0||jreference>=referenceReaderGdal.nrOfRow()){
785  if(rmse_opt[0]||regression_opt[0])
786  continue;
787  else{
788  cerr << jreference << " out of reference range!" << endl;
789  cerr << x << " " << y << " " << icol << " " << irow << endl;
790  cerr << x << " " << y << " " << ireference << " " << jreference << endl;
791  exit(1);
792  }
793  }
794  else{
795  referenceReaderGdal.readData(lineReference,static_cast<int>(jreference),band_opt[1]);
796  oldreferencerow=jreference;
797  }
798  }
799  bool flagged=false;
800  for(int iflag=0;iflag<nodata_opt.size();++iflag){
801  if((lineInput[icol]==nodata_opt[iflag])||(lineReference[ireference]==nodata_opt[iflag])){
802  if(output_opt.size())
803  lineOutput[icol]=nodata_opt[iflag];
804  flagged=true;
805  break;
806  }
807  }
808  if(mask_opt.size()){
809  maskReader.geo2image(x,y,imask,jmask);
810  if(jmask>=0&&jmask<maskReader.nrOfRow()){
811  if(jmask!=oldmaskrow)
812  maskReader.readData(lineMask,jmask);
813  for(int ivalue=0;ivalue<msknodata_opt.size();++ivalue){
814  if(lineMask[icol]==msknodata_opt[ivalue]){
815  flagged=true;
816  break;
817  }
818  }
819  }
820  }
821  if(!flagged){
822  if(rmse_opt[0]){//divide by image size to prevent overflow. At the end we need to take care about flagged pixels by normalizing...
823  rmse+=static_cast<double>(lineInput[icol]-lineReference[ireference])*(lineInput[icol]-lineReference[ireference])/inputReader.nrOfCol()/inputReader.nrOfRow();
824  }
825  else if(regression_opt[0]){
826  bufferInput.push_back(lineInput[icol]);
827  bufferReference.push_back(lineReference[ireference]);
828  }
829 
830  if(confusion_opt[0]){
831  ++ntotalValidation;
832  int rc=distance(referenceRange.begin(),find(referenceRange.begin(),referenceRange.end(),lineReference[ireference]));
833  int ic=distance(inputRange.begin(),find(inputRange.begin(),inputRange.end(),lineInput[icol]));
834  assert(rc<nclass);
835  assert(ic<nclass);
836  ++nvalidation[rc];
837  ++resultClass[rc][ic];
838  if(verbose_opt[0]>1)
839  cout << "increment: " << rc << " " << referenceRange[rc] << " " << ic << " " << inputRange[ic] << endl;
840  cm.incrementResult(cm.getClass(rc),cm.getClass(ic),1);
841  }
842  if(lineInput[icol]==lineReference[ireference]){//correct
843  if(output_opt.size()){
844  lineOutput[icol]=valueE_opt[0];
845  if(nodata_opt.size()){
846  if(valueE_opt[0]==nodata_opt[0])
847  lineOutput[icol]=lineInput[icol];
848  }
849  }
850  }
851  else{//error
852  if(output_opt.empty()&&!confusion_opt[0]&&!rmse_opt[0]&&!regression_opt[0]){
853  isDifferent=true;
854  break;
855  }
856  if(output_opt.size()){
857  if(lineInput[icol]>lineReference[ireference])
858  lineOutput[icol]=valueO_opt[0];//omission error
859  else
860  lineOutput[icol]=valueC_opt[0];//commission error
861  }
862  }
863  }
864  else{
865  ++nflagged;
866  if(output_opt.size()){
867  if(nodata_opt.size())
868  lineOutput[icol]=nodata_opt[0];
869  else //should never occur?
870  lineOutput[icol]=0;
871  }
872  }
873  }
874  if(output_opt.size()){
875  try{
876  gdalWriter.writeData(lineOutput,irow);
877  }
878  catch(string errorstring){
879  cerr << "lineOutput.size(): " << lineOutput.size() << endl;
880  cerr << "gdalWriter.nrOfCol(): " << gdalWriter.nrOfCol() << endl;
881  cerr << errorstring << endl;
882  exit(1);
883  }
884  }
885  else if(isDifferent&&!confusion_opt[0]&&!rmse_opt[0]&&!regression_opt[0]){//we can break off here, files are different...
886  if(!verbose_opt[0])
887  pfnProgress(1.0,pszMessage,pProgressArg);
888  break;
889  }
890  progress=static_cast<float>(irow+1.0)/inputReader.nrOfRow();
891  if(!verbose_opt[0])
892  pfnProgress(progress,pszMessage,pProgressArg);
893  }
894  if(output_opt.size())
895  gdalWriter.close();
896  else if(!confusion_opt[0]){
897  if(rmse_opt[0]){
898  double normalization=1.0*inputReader.nrOfCol()*inputReader.nrOfRow()/(inputReader.nrOfCol()*inputReader.nrOfRow()-nflagged);
899  if(verbose_opt[0]){
900  cout << "normalization: " << normalization << endl;
901  cout << "rmse before sqrt and normalization: " << rmse << endl;
902  }
903  cout << "--rmse " << sqrt(rmse/normalization) << endl;
904  }
905  else if(regression_opt[0]){
906  double err=0;
907  double c0=0;
908  double c1=1;
910  if(bufferInput.size()&&bufferReference.size()){
911  err=stat.linear_regression_err(bufferInput,bufferReference,c0,c1);
912  }
913  if(verbose_opt[0]){
914  cout << "bufferInput.size(): " << bufferInput.size() << endl;
915  cout << "bufferReference.size(): " << bufferReference.size() << endl;
916  double theMin=0;
917  double theMax=0;
918  stat.minmax(bufferInput,bufferInput.begin(),bufferInput.end(),theMin,theMax);
919  cout << "min, max input: " << theMin << ", " << theMax << endl;
920  theMin=0;
921  theMax=0;
922  stat.minmax(bufferReference,bufferReference.begin(),bufferReference.end(),theMin,theMax);
923  cout << "min, max reference: " << theMin << ", " << theMax << endl;
924  }
925  cout << "--c0 " << c0 << "--c1 " << c1 << " --rmse: " << err << endl;
926 
927  }
928  else if(isDifferent)
929  cout << input_opt[0] << " and " << reference_opt[0] << " are different" << endl;
930  else
931  cout << input_opt[0] << " and " << reference_opt[0] << " are identical" << endl;
932  }
933  referenceReaderGdal.close();
934  inputReader.close();
935  if(mask_opt.size())
936  maskReader.close();
937  }//raster dataset
938 
939  if(confusion_opt[0]){
940  cm.setFormat(cmformat_opt[0]);
941  cm.reportSE95(se95_opt[0]);
942  ofstream outputFile;
943  if(cmoutput_opt.size()){
944  outputFile.open(cmoutput_opt[0].c_str(),ios::out);
945  outputFile << cm << endl;
946  }
947  else
948  cout << cm << endl;
949  // cout << "class #samples userAcc prodAcc" << endl;
950  // double se95_ua=0;
951  // double se95_pa=0;
952  // double se95_oa=0;
953  // double dua=0;
954  // double dpa=0;
955  // double doa=0;
956  // for(int iclass=0;iclass<cm.nClasses();++iclass){
957  // dua=cm.ua_pct(classNames[iclass],&se95_ua);
958  // dpa=cm.pa_pct(classNames[iclass],&se95_pa);
959  // cout << cm.getClass(iclass) << " " << cm.nReference(cm.getClass(iclass)) << " " << dua << " (" << se95_ua << ")" << " " << dpa << " (" << se95_pa << ")" << endl;
960  // }
961  // doa=cm.oa(&se95_oa);
962  // cout << "Kappa: " << cm.kappa() << endl;
963  // cout << "Overall Accuracy: " << 100*doa << " (" << 100*se95_oa << ")" << endl;
964  }
965 }
STL namespace.
bool geo2image(double x, double y, double &i, double &j) const
Convert georeferenced coordinates (x and y) to image coordinates (column and row) ...
void close(void)
Set the memory (in MB) to cache a number of rows in memory.
int nrOfCol(void) const
Get the number of columns of this dataset.
Definition: ImgRasterGdal.h:98
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)
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
int nrOfRow(void) const
Get the number of rows of this dataset.
bool image2geo(double i, double j, double &x, double &y) const
Convert image coordinates (column and row) to georeferenced coordinates (x and y) ...
GDALColorTable * getColorTable(int band=0) const
Get the GDAL color table for this dataset as an instance of the GDALColorTable class.
void getRange(std::vector< short > &range, int Band=0)
Calculate the range of cell values in the image for a specific band (start counting from 0)...
bool isGeoRef() const
Is this dataset georeferenced (pixel size in y must be negative) ?
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...
std::string getProjection(void) const
Get the projection string (deprecated, use getProjectionRef instead)
GDALDataType getDataType(int band=0) const
Get the GDAL datatype for this dataset.
void close(void)
Close the image.
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.
std::string getInterleave() const
Get the band coding (interleave)