pktools 2.6.7
Processing Kernel for geospatial data
pkdiff.cc
1/**********************************************************************
2pkdiff.cc: program to compare two raster image files
3Copyright (C) 2008-2014 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 <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/******************************************************************************/
93using namespace std;
94
95int 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}
int nrOfRow(void) const
Get the number of rows of this dataset.
bool isGeoRef() const
Is this dataset georeferenced (pixel size in y must be negative) ?
bool geo2image(double x, double y, double &i, double &j) const
Convert georeferenced coordinates (x and y) to image coordinates (column and row)
std::string getProjection(void) const
Get the projection string (deprecated, use getProjectionRef instead)
GDALColorTable * getColorTable(int band=0) const
Get the GDAL color table for this dataset as an instance of the GDALColorTable class.
void copyGeoTransform(const ImgRasterGdal &imgSrc)
Copy geotransform information from another georeferenced image.
int nrOfCol(void) const
Get the number of columns of this dataset.
Definition: ImgRasterGdal.h:98
std::string getInterleave() const
Get the band coding (interleave)
bool image2geo(double i, double j, double &x, double &y) const
Convert image coordinates (column and row) to georeferenced coordinates (x and y)
CPLErr GDALSetNoDataValue(double noDataValue, int band=0)
Set the GDAL (internal) no data value for this data set. Only a single no data value per band is supp...
GDALDataType getDataType(int band=0) const
Get the GDAL datatype for this dataset.
void readData(T &value, int col, int row, int band=0)
Read a single pixel cell value at a specific column and row for a specific band (all indices start co...
Definition: ImgReaderGdal.h:95
void close(void)
Set the memory (in MB) to cache a number of rows in memory.
void open(const std::string &filename, const GDALAccess &readMode=GA_ReadOnly)
Open an image.
void 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).
void close(void)
Close the image.
void open(const std::string &filename, const ImgReaderGdal &imgSrc, const std::vector< std::string > &options=std::vector< std::string >())
Open an image for writing, copying image attributes from a source image. Image is directly written to...
void setColorTable(const std::string &filename, int band=0)
Set the color table using an (ASCII) file with 5 columns (value R G B alpha)
bool writeData(T &value, int col, int row, int band=0)
Write a single pixel cell value at a specific column and row for a specific band (all indices start c...
Definition: ImgWriterGdal.h:96