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" 95 int main(
int argc,
char *argv[])
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");
112 Optionpk<short> classvalue_opt(
"r",
"reclass",
"List of class values (use same order as in classname option).");
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);
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);
134 colorTable_opt.setHide(1);
135 option_opt.setHide(1);
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);
164 valueE_opt.retrieveOption(argc,argv);
165 valueO_opt.retrieveOption(argc,argv);
166 valueC_opt.retrieveOption(argc,argv);
167 verbose_opt.retrieveOption(argc,argv);
169 catch(
string predefinedString){
170 std::cout << predefinedString << std::endl;
175 cout <<
"Usage: pkdiff -i input -ref reference" << endl;
177 std::cout <<
"short option -h shows basic options only, use long option --help to show all options" << std::endl;
185 cout <<
"flag(s) set to";
186 for(
int iflag=0;iflag<nodata_opt.size();++iflag)
187 cout <<
" " << nodata_opt[iflag];
191 if(input_opt.empty()){
192 std::cerr <<
"No input file provided (use option -i). Use --help for help information" << std::endl;
195 if(reference_opt.empty()){
196 std::cerr <<
"No reference file provided (use option -ref). Use --help for help information" << std::endl;
202 if(band_opt.size()<2)
203 band_opt.push_back(band_opt[0]);
206 while(mask_opt.size()<input_opt.size())
207 mask_opt.push_back(mask_opt[0]);
208 vector<short> inputRange;
209 vector<short> referenceRange;
212 map<string,short> classValueMap;
213 vector<std::string> nameVector(255);
214 vector<string> classNames;
216 unsigned int ntotalValidation=0;
217 unsigned int nflagged=0;
220 vector<float> producer;
221 vector<unsigned int> nvalidation;
223 if(confusion_opt[0]){
231 cout <<
"opening input image file " << input_opt[0] << endl;
232 inputReader.
open(input_opt[0]);
235 cerr << error << endl;
238 inputReader.
getRange(inputRange,band_opt[0]);
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);
248 nclass=inputRange.size();
250 cout <<
"nclass (inputRange.size()): " << nclass << endl;
251 cout <<
"input range: " << endl;
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];
262 for(
int rc=0;rc<inputRange.size();++rc){
263 classNames.push_back(type2string(inputRange[rc]));
265 cout << inputRange[rc] << endl;
267 cm.setClassNames(classNames);
269 cout <<
"class names: " << endl;
270 for(
int iclass=0;iclass<cm.nClasses();++iclass)
271 cout << iclass <<
" " << cm.getClass(iclass) << endl;
273 resultClass.resize(nclass,nclass);
275 producer.resize(nclass);
276 nvalidation.resize(nclass);
278 for(
int rc=0;rc<nclass;++rc){
279 for(
int ic=0;ic<nclass;++ic)
280 resultClass[rc][ic]=0;
285 bool isDifferent=
false;
286 bool refIsRaster=
false;
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();
296 std::cout <<
"reference is vector dataset" << std::endl;
298 catch(
string errorString){
302 std::cout <<
"reference is raster dataset" << std::endl;
304 const char* pszMessage;
305 void* pProgressArg=NULL;
306 GDALProgressFunc pfnProgress=GDALTermProgress;
310 for(
int iinput=0;iinput<input_opt.size();++iinput){
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;
319 inputReader.
open(input_opt[iinput]);
321 maskReader.
open(mask_opt[iinput]);
325 referenceReaderOgr.open(reference_opt[iref]);
328 cerr << error << endl;
332 referenceRange=inputRange;
335 if(output_opt.size()){
337 ogrWriter.open(output_opt[iref],ogrformat_opt[0]);
340 cerr << error << endl;
344 int nlayer=referenceReaderOgr.getDataSource()->GetLayerCount();
345 for(
int ilayer=0;ilayer<nlayer;++ilayer){
347 OGRLayer *readLayer=referenceReaderOgr.getLayer(ilayer);
349 string currentLayername=readLayer->GetName();
351 if(find(layer_opt.begin(),layer_opt.end(),currentLayername)==layer_opt.end())
354 pfnProgress(progress,pszMessage,pProgressArg);
356 cout <<
"processing layer " << readLayer->GetName() << endl;
358 readLayer->ResetReading();
359 OGRLayer *writeLayer;
360 if(output_opt.size()){
362 cout <<
"creating output vector file " << output_opt[0] << endl;
364 char **papszOptions=NULL;
366 cout <<
"creating layer: " << readLayer->GetName() << endl;
368 writeLayer=ogrWriter.createLayer(readLayer->GetName(), referenceReaderOgr.getProjection(ilayer), wkbPoint, papszOptions);
371 cout <<
"created layer" << endl;
372 cout <<
"copy fields from " << reference_opt[iref] << endl;
374 ogrWriter.copyFields(referenceReaderOgr,ilayer,ilayer);
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)))
383 fs << labelclass_opt[0] <<
"_" << windowJ <<
"_" << windowI;
385 fs << labelclass_opt[0];
387 cout <<
"creating field " << fs.str() << endl;
388 ogrWriter.createField(fs.str(),OFTInteger,ilayer);
392 OGRFeature *readFeature;
393 OGRFeature *writeFeature;
395 unsigned int nfeatureInLayer=readLayer->GetFeatureCount();
396 unsigned int ifeature=0;
397 while( (readFeature = readLayer->GetNextFeature()) != NULL ){
399 cout <<
"sample " << ++isample << endl;
402 OGRGeometry *poGeometry;
403 OGRPoint centroidPoint;
405 poGeometry = readFeature->GetGeometryRef();
409 else if(wkbFlatten(poGeometry->getGeometryType()) == wkbMultiPolygon){
410 OGRMultiPolygon readPolygon = *((OGRMultiPolygon *) poGeometry);
411 readPolygon = *((OGRMultiPolygon *) poGeometry);
412 readPolygon.Centroid(¢roidPoint);
413 poPoint=¢roidPoint;
415 else if(wkbFlatten(poGeometry->getGeometryType()) == wkbPolygon){
416 OGRPolygon readPolygon=*((OGRPolygon *) poGeometry);
417 readPolygon.Centroid(¢roidPoint);
418 poPoint=¢roidPoint;
420 else if(wkbFlatten(poGeometry->getGeometryType()) == wkbPoint )
421 poPoint = (OGRPoint *) poGeometry;
423 std::cerr <<
"Warning: skipping feature (not of type point or polygon)" << std::endl;
429 vector<double> inputValues;
430 bool isHomogeneous=
true;
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];
441 referenceValue=readFeature->GetFieldAsInteger(readFeature->GetFieldIndex(labelref_opt[0].c_str()));
443 cout <<
"reference value: " << referenceValue << endl;
445 bool pixelFlagged=
false;
446 bool maskFlagged=
false;
447 for(
int iflag=0;iflag<nodata_opt.size();++iflag){
448 if(referenceValue==nodata_opt[iflag])
453 double i_centre,j_centre;
455 inputReader.
geo2image(x,y,i_centre,j_centre);
461 j_centre=
static_cast<int>(j_centre);
462 i_centre=
static_cast<int>(i_centre);
464 if(static_cast<int>(j_centre)<0||
static_cast<int>(j_centre)>=inputReader.
nrOfRow())
467 if(static_cast<int>(i_centre)<0||
static_cast<int>(i_centre)>=inputReader.
nrOfCol())
470 if(output_opt.size()){
471 writeFeature = OGRFeature::CreateFeature(writeLayer->GetLayerDefn());
473 int nfield=readFeature->GetFieldCount();
474 writeFeature->SetGeometry(poPoint);
476 cout <<
"copying fields from " << reference_opt[0] << endl;
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]));
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)))
496 int j=j_centre+windowJ;
498 if(static_cast<int>(j)<0||
static_cast<int>(j)>=inputReader.
nrOfRow())
500 int i=i_centre+windowI;
502 if(static_cast<int>(i)<0||
static_cast<int>(i)>=inputReader.
nrOfCol())
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()))
511 cout <<
"input value: " << inputValue << endl;
513 for(
int iflag=0;iflag<nodata_opt.size();++iflag){
514 if(inputValue==nodata_opt[iflag]){
521 maskReader.
readData(maskValue,i,j,0);
522 for(
int ivalue=0;ivalue<msknodata_opt.size();++ivalue){
523 if(msknodata_opt[ivalue]>=0){
524 if(maskValue==msknodata_opt[ivalue]){
530 if(maskValue!=-msknodata_opt[ivalue])
539 pixelFlagged=pixelFlagged||maskFlagged;
543 windowAllFlagged=
false;
548 if(homogeneous_opt[0]){
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]){
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);
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)));
569 ++resultClass[rc][ic];
571 cout <<
"increment: " << rc <<
" " << referenceRange[rc] <<
" " << ic <<
" " << inputRange[ic] << endl;
572 cm.incrementResult(cm.getClass(rc),cm.getClass(ic),1);
575 if(inputValue==referenceValue){
576 outputValue=valueE_opt[0];
577 if(nodata_opt.size()){
578 if(valueE_opt[0]==nodata_opt[0])
579 outputValue=inputValue;
582 else if(inputValue>referenceValue)
583 outputValue=valueO_opt[0];
585 outputValue=valueC_opt[0];
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)))
593 int j=j_centre+windowJ;
595 if(static_cast<int>(j)<0||
static_cast<int>(j)>=inputReader.
nrOfRow())
597 int i=i_centre+windowI;
599 if(static_cast<int>(i)<0||
static_cast<int>(i)>=inputReader.
nrOfCol())
601 if(!windowAllFlagged){
604 fs << labelclass_opt[0] <<
"_" << windowJ <<
"_" << windowI;
606 fs << labelclass_opt[0];
607 if(output_opt.size())
608 writeFeature->SetField(fs.str().c_str(),inputValue);
609 if(!windowJ&&!windowI){
610 if(confusion_opt[0]){
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);
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)));
627 ++resultClass[rc][ic];
629 cout <<
"increment: " << rc <<
" " << referenceRange[rc] <<
" " << ic <<
" " << inputRange[ic] << endl;
630 cm.incrementResult(cm.getClass(rc),cm.getClass(ic),1);
633 if(inputValue==referenceValue){
634 outputValue=valueE_opt[0];
635 if(nodata_opt.size()){
636 if(valueE_opt[0]==nodata_opt[0])
637 outputValue=inputValue;
640 else if(inputValue>referenceValue)
641 outputValue=valueO_opt[0];
643 outputValue=valueC_opt[0];
649 if(output_opt.size()){
650 if(!windowAllFlagged){
652 cout <<
"creating feature" << endl;
653 if(writeLayer->CreateFeature( writeFeature ) != OGRERR_NONE ){
654 string errorString=
"Failed to create feature in OGR vector file";
658 OGRFeature::DestroyFeature( writeFeature );
661 progress=
static_cast<float>(ifeature+1)/nfeatureInLayer;
662 pfnProgress(progress,pszMessage,pProgressArg);
665 if(output_opt.size())
667 referenceReaderOgr.close();
673 pfnProgress(1.0,pszMessage,pProgressArg);
678 inputReader.
open(input_opt[0]);
680 maskReader.
open(mask_opt[0]);
681 if(output_opt.size()){
683 cout <<
"opening output image " << output_opt[0] << endl;
684 if(option_opt.findSubstring(
"INTERLEAVE=")==option_opt.end()){
685 string theInterleave=
"INTERLEAVE=";
687 option_opt.push_back(theInterleave);
690 if(nodata_opt.size())
693 if(colorTable_opt.size())
697 cout <<
"set colortable from input image" << endl;
701 else if(verbose_opt[0])
702 cout <<
"no output image defined" << endl;
706 cout << error << endl;
710 vector<double> lineInput(inputReader.
nrOfCol());
711 vector<double> lineMask(maskReader.
nrOfCol());
712 vector<double> lineOutput;
713 vector<double> bufferInput;
714 vector<double> bufferReference;
715 if(output_opt.size())
716 lineOutput.resize(inputReader.
nrOfCol());
720 double oldreferencerow=-1;
721 double oldmaskrow=-1;
724 referenceReaderGdal.
open(reference_opt[0]);
727 cerr << error << endl;
731 assert(referenceReaderGdal.
isGeoRef());
733 cerr <<
"Warning: projection of input image and reference image are different" << endl;
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);
745 cout <<
"reference range: " << endl;
746 for(
int rc=0;rc<referenceRange.size();++rc)
747 cout << referenceRange[rc] << endl;
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;
756 cout <<
"reference range is not equal to input range!" << endl;
757 cout << input_opt[0] <<
" and " << reference_opt[0] <<
" are different" << endl;
763 for(irow=0;irow<inputReader.
nrOfRow();++irow){
765 inputReader.
readData(lineInput,irow,band_opt[0]);
767 double ireference,jreference;
769 for(icol=0;icol<inputReader.
nrOfCol();++icol){
772 referenceReaderGdal.
geo2image(x,y,ireference,jreference);
773 if(ireference<0||ireference>=referenceReaderGdal.
nrOfCol()){
774 if(rmse_opt[0]||regression_opt[0])
777 cerr << ireference <<
" out of reference range!" << endl;
778 cerr << x <<
" " << y <<
" " << icol <<
" " << irow << endl;
779 cerr << x <<
" " << y <<
" " << ireference <<
" " << jreference << endl;
783 if(jreference!=oldreferencerow){
784 if(jreference<0||jreference>=referenceReaderGdal.
nrOfRow()){
785 if(rmse_opt[0]||regression_opt[0])
788 cerr << jreference <<
" out of reference range!" << endl;
789 cerr << x <<
" " << y <<
" " << icol <<
" " << irow << endl;
790 cerr << x <<
" " << y <<
" " << ireference <<
" " << jreference << endl;
795 referenceReaderGdal.
readData(lineReference,static_cast<int>(jreference),band_opt[1]);
796 oldreferencerow=jreference;
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];
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]){
823 rmse+=
static_cast<double>(lineInput[icol]-lineReference[ireference])*(lineInput[icol]-lineReference[ireference])/inputReader.
nrOfCol()/inputReader.
nrOfRow();
825 else if(regression_opt[0]){
826 bufferInput.push_back(lineInput[icol]);
827 bufferReference.push_back(lineReference[ireference]);
830 if(confusion_opt[0]){
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]));
837 ++resultClass[rc][ic];
839 cout <<
"increment: " << rc <<
" " << referenceRange[rc] <<
" " << ic <<
" " << inputRange[ic] << endl;
840 cm.incrementResult(cm.getClass(rc),cm.getClass(ic),1);
842 if(lineInput[icol]==lineReference[ireference]){
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];
852 if(output_opt.empty()&&!confusion_opt[0]&&!rmse_opt[0]&&!regression_opt[0]){
856 if(output_opt.size()){
857 if(lineInput[icol]>lineReference[ireference])
858 lineOutput[icol]=valueO_opt[0];
860 lineOutput[icol]=valueC_opt[0];
866 if(output_opt.size()){
867 if(nodata_opt.size())
868 lineOutput[icol]=nodata_opt[0];
874 if(output_opt.size()){
878 catch(
string errorstring){
879 cerr <<
"lineOutput.size(): " << lineOutput.size() << endl;
880 cerr <<
"gdalWriter.nrOfCol(): " << gdalWriter.
nrOfCol() << endl;
881 cerr << errorstring << endl;
885 else if(isDifferent&&!confusion_opt[0]&&!rmse_opt[0]&&!regression_opt[0]){
887 pfnProgress(1.0,pszMessage,pProgressArg);
890 progress=
static_cast<float>(irow+1.0)/inputReader.
nrOfRow();
892 pfnProgress(progress,pszMessage,pProgressArg);
894 if(output_opt.size())
896 else if(!confusion_opt[0]){
900 cout <<
"normalization: " << normalization << endl;
901 cout <<
"rmse before sqrt and normalization: " << rmse << endl;
903 cout <<
"--rmse " << sqrt(rmse/normalization) << endl;
905 else if(regression_opt[0]){
910 if(bufferInput.size()&&bufferReference.size()){
911 err=stat.linear_regression_err(bufferInput,bufferReference,c0,c1);
914 cout <<
"bufferInput.size(): " << bufferInput.size() << endl;
915 cout <<
"bufferReference.size(): " << bufferReference.size() << endl;
918 stat.minmax(bufferInput,bufferInput.begin(),bufferInput.end(),theMin,theMax);
919 cout <<
"min, max input: " << theMin <<
", " << theMax << endl;
922 stat.minmax(bufferReference,bufferReference.begin(),bufferReference.end(),theMin,theMax);
923 cout <<
"min, max reference: " << theMin <<
", " << theMax << endl;
925 cout <<
"--c0 " << c0 <<
"--c1 " << c1 <<
" --rmse: " << err << endl;
929 cout << input_opt[0] <<
" and " << reference_opt[0] <<
" are different" << endl;
931 cout << input_opt[0] <<
" and " << reference_opt[0] <<
" are identical" << endl;
933 referenceReaderGdal.
close();
939 if(confusion_opt[0]){
940 cm.setFormat(cmformat_opt[0]);
941 cm.reportSE95(se95_opt[0]);
943 if(cmoutput_opt.size()){
944 outputFile.open(cmoutput_opt[0].c_str(),ios::out);
945 outputFile << cm << endl;
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.
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...
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...
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)