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  referenceReaderOgr.close();
292  }
293  catch(string errorString){
294  //todo: sampleIsRaster will not work from GDAL 2.0!!?? (unification of driver for raster and vector datasets)
295  refIsRaster=true;
296  }
297  const char* pszMessage;
298  void* pProgressArg=NULL;
299  GDALProgressFunc pfnProgress=GDALTermProgress;
300  float progress=0;
301  // if(reference_opt[0].find(".shp")!=string::npos){
302  if(!refIsRaster){
303  for(int iinput=0;iinput<input_opt.size();++iinput){
304  if(verbose_opt[0])
305  cout << "Processing input " << input_opt[iinput] << endl;
306  if(output_opt.size())
307  assert(reference_opt.size()==output_opt.size());
308  for(int iref=0;iref<reference_opt.size();++iref){
309  cout << "reference " << reference_opt[iref] << endl;
310  // assert(reference_opt[iref].find(".shp")!=string::npos);
311  try{
312  inputReader.open(input_opt[iinput]);//,imagicX_opt[0],imagicY_opt[0]);
313  if(mask_opt.size()){
314  maskReader.open(mask_opt[iinput]);
315  assert(inputReader.nrOfCol()==maskReader.nrOfCol());
316  assert(inputReader.nrOfRow()==maskReader.nrOfRow());
317  }
318  referenceReaderOgr.open(reference_opt[iref]);
319  }
320  catch(string error){
321  cerr << error << endl;
322  exit(1);
323  }
324  if(confusion_opt[0])
325  referenceRange=inputRange;
326 
327  ImgWriterOgr ogrWriter;
328  if(output_opt.size()){
329  try{
330  ogrWriter.open(output_opt[iref],ogrformat_opt[0]);
331  }
332  catch(string error){
333  cerr << error << endl;
334  exit(1);
335  }
336  }
337  int nlayer=referenceReaderOgr.getDataSource()->GetLayerCount();
338  for(int ilayer=0;ilayer<nlayer;++ilayer){
339  progress=0;
340  OGRLayer *readLayer=referenceReaderOgr.getLayer(ilayer);
341  // readLayer = referenceReaderOgr.getDataSource()->GetLayer(ilayer);
342  string currentLayername=readLayer->GetName();
343  if(layer_opt.size())
344  if(find(layer_opt.begin(),layer_opt.end(),currentLayername)==layer_opt.end())
345  continue;
346  if(!verbose_opt[0])
347  pfnProgress(progress,pszMessage,pProgressArg);
348  else
349  cout << "processing layer " << readLayer->GetName() << endl;
350 
351  readLayer->ResetReading();
352  OGRLayer *writeLayer;
353  if(output_opt.size()){
354  if(verbose_opt[0])
355  cout << "creating output vector file " << output_opt[0] << endl;
356  // assert(output_opt[0].find(".shp")!=string::npos);
357  char **papszOptions=NULL;
358  if(verbose_opt[0])
359  cout << "creating layer: " << readLayer->GetName() << endl;
360  // if(ogrWriter.createLayer(layername, referenceReaderOgr.getProjection(ilayer), referenceReaderOgr.getGeometryType(ilayer), papszOptions)==NULL)
361  writeLayer=ogrWriter.createLayer(readLayer->GetName(), referenceReaderOgr.getProjection(ilayer), wkbPoint, papszOptions);
362  assert(writeLayer);
363  if(verbose_opt[0]){
364  cout << "created layer" << endl;
365  cout << "copy fields from " << reference_opt[iref] << endl;
366  }
367  ogrWriter.copyFields(referenceReaderOgr,ilayer,ilayer);
368  //create extra field for classified label
369  short theDim=boundary_opt[0];
370  for(int windowJ=-theDim/2;windowJ<(theDim+1)/2;++windowJ){
371  for(int windowI=-theDim/2;windowI<(theDim+1)/2;++windowI){
372  if(disc_opt[0]&&(windowI*windowI+windowJ*windowJ>(theDim/2)*(theDim/2)))
373  continue;
374  ostringstream fs;
375  if(theDim>1)
376  fs << labelclass_opt[0] << "_" << windowJ << "_" << windowI;
377  else
378  fs << labelclass_opt[0];
379  if(verbose_opt[0])
380  cout << "creating field " << fs.str() << endl;
381  ogrWriter.createField(fs.str(),OFTInteger,ilayer);
382  }
383  }
384  }
385  OGRFeature *readFeature;
386  OGRFeature *writeFeature;
387  int isample=0;
388  unsigned int nfeatureInLayer=readLayer->GetFeatureCount();
389  unsigned int ifeature=0;
390  while( (readFeature = readLayer->GetNextFeature()) != NULL ){
391  if(verbose_opt[0])
392  cout << "sample " << ++isample << endl;
393  //get x and y from readFeature
394  double x,y;
395  OGRGeometry *poGeometry;
396  OGRPoint centroidPoint;
397  OGRPoint *poPoint;
398  poGeometry = readFeature->GetGeometryRef();
399  // assert( poGeometry != NULL && wkbFlatten(poGeometry->getGeometryType()) == wkbPoint );
400  if(poGeometry==NULL)
401  continue;
402  else if(wkbFlatten(poGeometry->getGeometryType()) == wkbMultiPolygon){
403  OGRMultiPolygon readPolygon = *((OGRMultiPolygon *) poGeometry);
404  readPolygon = *((OGRMultiPolygon *) poGeometry);
405  readPolygon.Centroid(&centroidPoint);
406  poPoint=&centroidPoint;
407  }
408  else if(wkbFlatten(poGeometry->getGeometryType()) == wkbPolygon){
409  OGRPolygon readPolygon=*((OGRPolygon *) poGeometry);
410  readPolygon.Centroid(&centroidPoint);
411  poPoint=&centroidPoint;
412  }
413  else if(wkbFlatten(poGeometry->getGeometryType()) == wkbPoint )
414  poPoint = (OGRPoint *) poGeometry;
415  else{
416  std::cerr << "Warning: skipping feature (not of type point or polygon)" << std::endl;
417  continue;
418  }
419  x=poPoint->getX();
420  y=poPoint->getY();
421  double inputValue;
422  vector<double> inputValues;
423  bool isHomogeneous=true;
424  short maskValue;
425  short outputValue;
426  //read referenceValue from feature
427  unsigned short referenceValue;
428  string referenceClassName;
429  if(classValueMap.size()){
430  referenceClassName=readFeature->GetFieldAsString(readFeature->GetFieldIndex(labelref_opt[0].c_str()));
431  referenceValue=classValueMap[referenceClassName];
432  }
433  else
434  referenceValue=readFeature->GetFieldAsInteger(readFeature->GetFieldIndex(labelref_opt[0].c_str()));
435  if(verbose_opt[0])
436  cout << "reference value: " << referenceValue << endl;
437 
438  bool pixelFlagged=false;
439  bool maskFlagged=false;
440  for(int iflag=0;iflag<nodata_opt.size();++iflag){
441  if(referenceValue==nodata_opt[iflag])
442  pixelFlagged=true;
443  }
444  if(pixelFlagged)
445  continue;
446  double i_centre,j_centre;
447  //input reader is georeferenced!
448  inputReader.geo2image(x,y,i_centre,j_centre);
449  // else{
450  // i_centre=x;
451  // j_centre=y;
452  // }
453  //nearest neighbour
454  j_centre=static_cast<int>(j_centre);
455  i_centre=static_cast<int>(i_centre);
456  //check if j_centre is out of bounds
457  if(static_cast<int>(j_centre)<0||static_cast<int>(j_centre)>=inputReader.nrOfRow())
458  continue;
459  //check if i_centre is out of bounds
460  if(static_cast<int>(i_centre)<0||static_cast<int>(i_centre)>=inputReader.nrOfCol())
461  continue;
462 
463  if(output_opt.size()){
464  writeFeature = OGRFeature::CreateFeature(writeLayer->GetLayerDefn());
465  assert(readFeature);
466  int nfield=readFeature->GetFieldCount();
467  writeFeature->SetGeometry(poPoint);
468  if(verbose_opt[0])
469  cout << "copying fields from " << reference_opt[0] << endl;
470  assert(readFeature);
471  assert(writeFeature);
472  vector<int> panMap(nfield);
473  vector<int>::iterator panit=panMap.begin();
474  for(int ifield=0;ifield<nfield;++ifield)
475  panMap[ifield]=ifield;
476  writeFeature->SetFieldsFrom(readFeature,&(panMap[0]));
477  // if(writeFeature->SetFrom(readFeature)!= OGRERR_NONE)
478  // cerr << "writing feature failed" << endl;
479  // if(verbose_opt[0])
480  // cout << "feature written" << endl;
481  }
482  bool windowAllFlagged=true;
483  bool windowHasFlag=false;
484  short theDim=boundary_opt[0];
485  for(int windowJ=-theDim/2;windowJ<(theDim+1)/2;++windowJ){
486  for(int windowI=-theDim/2;windowI<(theDim+1)/2;++windowI){
487  if(disc_opt[0]&&(windowI*windowI+windowJ*windowJ>(theDim/2)*(theDim/2)))
488  continue;
489  int j=j_centre+windowJ;
490  //check if j is out of bounds
491  if(static_cast<int>(j)<0||static_cast<int>(j)>=inputReader.nrOfRow())
492  continue;
493  int i=i_centre+windowI;
494  //check if i is out of bounds
495  if(static_cast<int>(i)<0||static_cast<int>(i)>=inputReader.nrOfCol())
496  continue;
497  if(verbose_opt[0])
498  cout << setprecision(12) << "reading image value at x,y " << x << "," << y << " (" << i << "," << j << "), ";
499  inputReader.readData(inputValue,i,j,band_opt[0]);
500  inputValues.push_back(inputValue);
501  if(inputValues.back()!=*(inputValues.begin()))
502  isHomogeneous=false;
503  if(verbose_opt[0])
504  cout << "input value: " << inputValue << endl;
505  pixelFlagged=false;
506  for(int iflag=0;iflag<nodata_opt.size();++iflag){
507  if(inputValue==nodata_opt[iflag]){
508  pixelFlagged=true;
509  break;
510  }
511  }
512  maskFlagged=false;//(msknodata_opt[ivalue]>=0)?false:true;
513  if(mask_opt.size()){
514  maskReader.readData(maskValue,i,j,0);
515  for(int ivalue=0;ivalue<msknodata_opt.size();++ivalue){
516  if(msknodata_opt[ivalue]>=0){//values set in msknodata_opt are invalid
517  if(maskValue==msknodata_opt[ivalue]){
518  maskFlagged=true;
519  break;
520  }
521  }
522  else{//only values set in msknodata_opt are valid
523  if(maskValue!=-msknodata_opt[ivalue])
524  maskFlagged=true;
525  else{
526  maskFlagged=false;
527  break;
528  }
529  }
530  }
531  }
532  pixelFlagged=pixelFlagged||maskFlagged;
533  if(pixelFlagged)
534  windowHasFlag=true;
535  else
536  windowAllFlagged=false;//at least one good pixel in neighborhood
537  }
538  }
539  //at this point we know the values for the entire window
540 
541  if(homogeneous_opt[0]){//only centre pixel
542  int j=j_centre;
543  int i=i_centre;
544  //flag if not all pixels are homogeneous or if at least one pixel flagged
545 
546  if(!windowHasFlag&&isHomogeneous){
547  if(output_opt.size())
548  writeFeature->SetField(labelclass_opt[0].c_str(),static_cast<int>(inputValue));
549  if(confusion_opt[0]){
550  ++ntotalValidation;
551  if(classValueMap.size()){
552  assert(inputValue<nameVector.size());
553  string className=nameVector[static_cast<unsigned short>(inputValue)];
554  cm.incrementResult(type2string<short>(classValueMap[referenceClassName]),type2string<short>(classValueMap[className]),1);
555  }
556  else{
557  int rc=distance(referenceRange.begin(),find(referenceRange.begin(),referenceRange.end(),static_cast<unsigned short>(referenceValue)));
558  int ic=distance(inputRange.begin(),find(inputRange.begin(),inputRange.end(),static_cast<unsigned short>(inputValue)));
559  assert(rc<nclass);
560  assert(ic<nclass);
561  ++nvalidation[rc];
562  ++resultClass[rc][ic];
563  if(verbose_opt[0]>1)
564  cout << "increment: " << rc << " " << referenceRange[rc] << " " << ic << " " << inputRange[ic] << endl;
565  cm.incrementResult(cm.getClass(rc),cm.getClass(ic),1);
566  }
567  }
568  if(inputValue==referenceValue){//correct
569  outputValue=valueE_opt[0];
570  if(nodata_opt.size()){
571  if(valueE_opt[0]==nodata_opt[0])
572  outputValue=inputValue;
573  }
574  }
575  else if(inputValue>referenceValue)//1=forest,2=non-forest
576  outputValue=valueO_opt[0];//omission error
577  else
578  outputValue=valueC_opt[0];//commission error
579  }
580  }
581  else{
582  for(int windowJ=-theDim/2;windowJ<(theDim+1)/2;++windowJ){
583  for(int windowI=-theDim/2;windowI<(theDim+1)/2;++windowI){
584  if(disc_opt[0]&&(windowI*windowI+windowJ*windowJ>(theDim/2)*(theDim/2)))
585  continue;
586  int j=j_centre+windowJ;
587  //check if j is out of bounds
588  if(static_cast<int>(j)<0||static_cast<int>(j)>=inputReader.nrOfRow())
589  continue;
590  int i=i_centre+windowI;
591  //check if i is out of bounds
592  if(static_cast<int>(i)<0||static_cast<int>(i)>=inputReader.nrOfCol())
593  continue;
594  if(!windowAllFlagged){
595  ostringstream fs;
596  if(theDim>1)
597  fs << labelclass_opt[0] << "_" << windowJ << "_" << windowI;
598  else
599  fs << labelclass_opt[0];
600  if(output_opt.size())
601  writeFeature->SetField(fs.str().c_str(),inputValue);
602  if(!windowJ&&!windowI){//centre pixel
603  if(confusion_opt[0]){
604  ++ntotalValidation;
605  if(classValueMap.size()){
606  assert(inputValue<nameVector.size());
607  string className=nameVector[static_cast<unsigned short>(inputValue)];
608  cm.incrementResult(type2string<short>(classValueMap[referenceClassName]),type2string<short>(classValueMap[className]),1);
609  }
610  else{
611  int rc=distance(referenceRange.begin(),find(referenceRange.begin(),referenceRange.end(),static_cast<unsigned short>(referenceValue)));
612  int ic=distance(inputRange.begin(),find(inputRange.begin(),inputRange.end(),static_cast<unsigned short>(inputValue)));
613  if(rc>=nclass)
614  continue;
615  if(ic>=nclass)
616  continue;
617  // assert(rc<nclass);
618  // assert(ic<nclass);
619  ++nvalidation[rc];
620  ++resultClass[rc][ic];
621  if(verbose_opt[0]>1)
622  cout << "increment: " << rc << " " << referenceRange[rc] << " " << ic << " " << inputRange[ic] << endl;
623  cm.incrementResult(cm.getClass(rc),cm.getClass(ic),1);
624  }
625  }
626  if(inputValue==referenceValue){//correct
627  outputValue=valueE_opt[0];
628  if(nodata_opt.size()){
629  if(valueE_opt[0]==nodata_opt[0])
630  outputValue=inputValue;
631  }
632  }
633  else if(inputValue>referenceValue)//1=forest,2=non-forest
634  outputValue=valueO_opt[0];//omission error
635  else
636  outputValue=valueC_opt[0];//commission error
637  }
638  }
639  }
640  }
641  }
642  if(output_opt.size()){
643  if(!windowAllFlagged){
644  if(verbose_opt[0])
645  cout << "creating feature" << endl;
646  if(writeLayer->CreateFeature( writeFeature ) != OGRERR_NONE ){
647  string errorString="Failed to create feature in OGR vector file";
648  throw(errorString);
649  }
650  }
651  OGRFeature::DestroyFeature( writeFeature );
652  }
653  ++ifeature;
654  progress=static_cast<float>(ifeature+1)/nfeatureInLayer;
655  pfnProgress(progress,pszMessage,pProgressArg);
656  }//next feature
657  }//next layer
658  if(output_opt.size())
659  ogrWriter.close();
660  referenceReaderOgr.close();
661  inputReader.close();
662  if(mask_opt.size())
663  maskReader.close();
664  }//next reference
665  }//next input
666  pfnProgress(1.0,pszMessage,pProgressArg);
667  }//reference is OGR vector
668  else{//reference is GDAL raster
669  ImgWriterGdal gdalWriter;
670  try{
671  inputReader.open(input_opt[0]);
672  if(mask_opt.size())
673  maskReader.open(mask_opt[0]);
674  if(output_opt.size()){
675  if(verbose_opt[0])
676  cout << "opening output image " << output_opt[0] << endl;
677  if(option_opt.findSubstring("INTERLEAVE=")==option_opt.end()){
678  string theInterleave="INTERLEAVE=";
679  theInterleave+=inputReader.getInterleave();
680  option_opt.push_back(theInterleave);
681  }
682  gdalWriter.open(output_opt[0],inputReader.nrOfCol(),inputReader.nrOfRow(),1,inputReader.getDataType(),oformat_opt[0],option_opt);
683  if(nodata_opt.size())
684  gdalWriter.GDALSetNoDataValue(nodata_opt[0]);
685  gdalWriter.copyGeoTransform(inputReader);
686  if(colorTable_opt.size())
687  gdalWriter.setColorTable(colorTable_opt[0]);
688  else if(inputReader.getColorTable()!=NULL){
689  if(verbose_opt[0])
690  cout << "set colortable from input image" << endl;
691  gdalWriter.setColorTable(inputReader.getColorTable());
692  }
693  }
694  else if(verbose_opt[0])
695  cout << "no output image defined" << endl;
696 
697  }
698  catch(string error){
699  cout << error << endl;
700  exit(2);
701  }
702  //todo: support different data types!
703  vector<double> lineInput(inputReader.nrOfCol());
704  vector<double> lineMask(maskReader.nrOfCol());
705  vector<double> lineOutput;
706  vector<double> bufferInput;//for regression
707  vector<double> bufferReference;//for regression
708  if(output_opt.size())
709  lineOutput.resize(inputReader.nrOfCol());
710 
711  int irow=0;
712  int icol=0;
713  double oldreferencerow=-1;
714  double oldmaskrow=-1;
715  ImgReaderGdal referenceReaderGdal;
716  try{
717  referenceReaderGdal.open(reference_opt[0]);//,rmagicX_opt[0],rmagicY_opt[0]);
718  }
719  catch(string error){
720  cerr << error << endl;
721  exit(1);
722  }
723  if(inputReader.isGeoRef()){
724  assert(referenceReaderGdal.isGeoRef());
725  if(inputReader.getProjection()!=referenceReaderGdal.getProjection())
726  cerr << "Warning: projection of input image and reference image are different" << endl;
727  }
728  vector<double> lineReference(referenceReaderGdal.nrOfCol());
729  if(confusion_opt[0]){
730  referenceReaderGdal.getRange(referenceRange,band_opt[1]);
731  for(int iflag=0;iflag<nodata_opt.size();++iflag){
732  vector<short>::iterator fit;
733  fit=find(referenceRange.begin(),referenceRange.end(),static_cast<unsigned short>(nodata_opt[iflag]));
734  if(fit!=referenceRange.end())
735  referenceRange.erase(fit);
736  }
737  if(verbose_opt[0]){
738  cout << "reference range: " << endl;
739  for(int rc=0;rc<referenceRange.size();++rc)
740  cout << referenceRange[rc] << endl;
741  }
742  if(referenceRange.size()!=inputRange.size()){
743  if(confusion_opt[0]||output_opt.size()){
744  cout << "reference range is not equal to input range!" << endl;
745  cout << "Kappa: " << 0 << endl;
746  cout << "total weighted: " << 0 << endl;
747  }
748  else
749  cout << "reference range is not equal to input range!" << endl;
750  cout << input_opt[0] << " and " << reference_opt[0] << " are different" << endl;
751  exit(1);
752  }
753  }
754  double rmse=0;
755  // for(irow=0;irow<inputReader.nrOfRow()&&!isDifferent;++irow){
756  for(irow=0;irow<inputReader.nrOfRow();++irow){
757  //read line in lineInput, lineReference and lineMask
758  inputReader.readData(lineInput,irow,band_opt[0]);
759  double x,y;//geo coordinates
760  double ireference,jreference;//image coordinates in reference image
761  double imask,jmask;//image coordinates in mask image
762  for(icol=0;icol<inputReader.nrOfCol();++icol){
763  //find col in reference
764  inputReader.image2geo(icol,irow,x,y);
765  referenceReaderGdal.geo2image(x,y,ireference,jreference);
766  if(ireference<0||ireference>=referenceReaderGdal.nrOfCol()){
767  if(rmse_opt[0]||regression_opt[0])
768  continue;
769  else{
770  cerr << ireference << " out of reference range!" << endl;
771  cerr << x << " " << y << " " << icol << " " << irow << endl;
772  cerr << x << " " << y << " " << ireference << " " << jreference << endl;
773  exit(1);
774  }
775  }
776  if(jreference!=oldreferencerow){
777  if(jreference<0||jreference>=referenceReaderGdal.nrOfRow()){
778  if(rmse_opt[0]||regression_opt[0])
779  continue;
780  else{
781  cerr << jreference << " out of reference range!" << endl;
782  cerr << x << " " << y << " " << icol << " " << irow << endl;
783  cerr << x << " " << y << " " << ireference << " " << jreference << endl;
784  exit(1);
785  }
786  }
787  else{
788  referenceReaderGdal.readData(lineReference,static_cast<int>(jreference),band_opt[1]);
789  oldreferencerow=jreference;
790  }
791  }
792  bool flagged=false;
793  for(int iflag=0;iflag<nodata_opt.size();++iflag){
794  if((lineInput[icol]==nodata_opt[iflag])||(lineReference[ireference]==nodata_opt[iflag])){
795  if(output_opt.size())
796  lineOutput[icol]=nodata_opt[iflag];
797  flagged=true;
798  break;
799  }
800  }
801  if(mask_opt.size()){
802  maskReader.geo2image(x,y,imask,jmask);
803  if(jmask>=0&&jmask<maskReader.nrOfRow()){
804  if(jmask!=oldmaskrow)
805  maskReader.readData(lineMask,jmask);
806  for(int ivalue=0;ivalue<msknodata_opt.size();++ivalue){
807  if(lineMask[icol]==msknodata_opt[ivalue]){
808  flagged=true;
809  break;
810  }
811  }
812  }
813  }
814  if(!flagged){
815  if(rmse_opt[0]){//divide by image size to prevent overflow. At the end we need to take care about flagged pixels by normalizing...
816  rmse+=static_cast<double>(lineInput[icol]-lineReference[ireference])*(lineInput[icol]-lineReference[ireference])/inputReader.nrOfCol()/inputReader.nrOfRow();
817  }
818  else if(regression_opt[0]){
819  bufferInput.push_back(lineInput[icol]);
820  bufferReference.push_back(lineReference[ireference]);
821  }
822 
823  if(confusion_opt[0]){
824  ++ntotalValidation;
825  int rc=distance(referenceRange.begin(),find(referenceRange.begin(),referenceRange.end(),lineReference[ireference]));
826  int ic=distance(inputRange.begin(),find(inputRange.begin(),inputRange.end(),lineInput[icol]));
827  assert(rc<nclass);
828  assert(ic<nclass);
829  ++nvalidation[rc];
830  ++resultClass[rc][ic];
831  if(verbose_opt[0]>1)
832  cout << "increment: " << rc << " " << referenceRange[rc] << " " << ic << " " << inputRange[ic] << endl;
833  cm.incrementResult(cm.getClass(rc),cm.getClass(ic),1);
834  }
835  if(lineInput[icol]==lineReference[ireference]){//correct
836  if(output_opt.size()){
837  lineOutput[icol]=valueE_opt[0];
838  if(nodata_opt.size()){
839  if(valueE_opt[0]==nodata_opt[0])
840  lineOutput[icol]=lineInput[icol];
841  }
842  }
843  }
844  else{//error
845  if(output_opt.empty()&&!confusion_opt[0]&&!rmse_opt[0]&&!regression_opt[0]){
846  isDifferent=true;
847  break;
848  }
849  if(output_opt.size()){
850  if(lineInput[icol]>lineReference[ireference])
851  lineOutput[icol]=valueO_opt[0];//omission error
852  else
853  lineOutput[icol]=valueC_opt[0];//commission error
854  }
855  }
856  }
857  else{
858  ++nflagged;
859  if(output_opt.size()){
860  if(nodata_opt.size())
861  lineOutput[icol]=nodata_opt[0];
862  else //should never occur?
863  lineOutput[icol]=0;
864  }
865  }
866  }
867  if(output_opt.size()){
868  try{
869  gdalWriter.writeData(lineOutput,irow);
870  }
871  catch(string errorstring){
872  cerr << "lineOutput.size(): " << lineOutput.size() << endl;
873  cerr << "gdalWriter.nrOfCol(): " << gdalWriter.nrOfCol() << endl;
874  cerr << errorstring << endl;
875  exit(1);
876  }
877  }
878  else if(isDifferent&&!confusion_opt[0]&&!rmse_opt[0]&&!regression_opt[0]){//we can break off here, files are different...
879  if(!verbose_opt[0])
880  pfnProgress(1.0,pszMessage,pProgressArg);
881  break;
882  }
883  progress=static_cast<float>(irow+1.0)/inputReader.nrOfRow();
884  if(!verbose_opt[0])
885  pfnProgress(progress,pszMessage,pProgressArg);
886  }
887  if(output_opt.size())
888  gdalWriter.close();
889  else if(!confusion_opt[0]){
890  if(rmse_opt[0]){
891  double normalization=1.0*inputReader.nrOfCol()*inputReader.nrOfRow()/(inputReader.nrOfCol()*inputReader.nrOfRow()-nflagged);
892  if(verbose_opt[0]){
893  cout << "normalization: " << normalization << endl;
894  cout << "rmse before sqrt and normalization: " << rmse << endl;
895  }
896  cout << "--rmse " << sqrt(rmse/normalization) << endl;
897  }
898  else if(regression_opt[0]){
899  double err=0;
900  double c0=0;
901  double c1=1;
903  if(bufferInput.size()&&bufferReference.size()){
904  err=stat.linear_regression_err(bufferInput,bufferReference,c0,c1);
905  }
906  if(verbose_opt[0]){
907  cout << "bufferInput.size(): " << bufferInput.size() << endl;
908  cout << "bufferReference.size(): " << bufferReference.size() << endl;
909  double theMin=0;
910  double theMax=0;
911  stat.minmax(bufferInput,bufferInput.begin(),bufferInput.end(),theMin,theMax);
912  cout << "min, max input: " << theMin << ", " << theMax << endl;
913  theMin=0;
914  theMax=0;
915  stat.minmax(bufferReference,bufferReference.begin(),bufferReference.end(),theMin,theMax);
916  cout << "min, max reference: " << theMin << ", " << theMax << endl;
917  }
918  cout << "--c0 " << c0 << "--c1 " << c1 << " --rmse: " << err << endl;
919 
920  }
921  else if(isDifferent)
922  cout << input_opt[0] << " and " << reference_opt[0] << " are different" << endl;
923  else
924  cout << input_opt[0] << " and " << reference_opt[0] << " are identical" << endl;
925  }
926  referenceReaderGdal.close();
927  inputReader.close();
928  if(mask_opt.size())
929  maskReader.close();
930  }//raster dataset
931 
932  if(confusion_opt[0]){
933  cm.setFormat(cmformat_opt[0]);
934  cm.reportSE95(se95_opt[0]);
935  ofstream outputFile;
936  if(cmoutput_opt.size()){
937  outputFile.open(cmoutput_opt[0].c_str(),ios::out);
938  outputFile << cm << endl;
939  }
940  else
941  cout << cm << endl;
942  // cout << "class #samples userAcc prodAcc" << endl;
943  // double se95_ua=0;
944  // double se95_pa=0;
945  // double se95_oa=0;
946  // double dua=0;
947  // double dpa=0;
948  // double doa=0;
949  // for(int iclass=0;iclass<cm.nClasses();++iclass){
950  // dua=cm.ua_pct(classNames[iclass],&se95_ua);
951  // dpa=cm.pa_pct(classNames[iclass],&se95_pa);
952  // cout << cm.getClass(iclass) << " " << cm.nReference(cm.getClass(iclass)) << " " << dua << " (" << se95_ua << ")" << " " << dpa << " (" << se95_pa << ")" << endl;
953  // }
954  // doa=cm.oa(&se95_oa);
955  // cout << "Kappa: " << cm.kappa() << endl;
956  // cout << "Overall Accuracy: " << 100*doa << " (" << 100*se95_oa << ")" << endl;
957  }
958 }
bool geo2image(double x, double y, double &i, double &j) const
Convert georeferenced coordinates (x and y) to image coordinates (column and row) ...
GDALColorTable * getColorTable(int band=0) const
Get the GDAL color table for this dataset as an instance of the GDALColorTable class.
STL namespace.
GDALDataType getDataType(int band=0) const
Get the GDAL datatype for this dataset.
bool isGeoRef() const
Is this dataset georeferenced (pixel size in y must be negative) ?
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
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 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 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