This is an automated email from the git hooks/post-receive script. sebastic pushed a commit to branch master in repository pktools.
commit f723b1d0f6343c8b4431d403b94896212217f2dc Author: Bas Couwenberg <[email protected]> Date: Sun Mar 8 11:59:39 2015 +0100 Imported Upstream version 2.6.3 --- ChangeLog | 16 + configure | 32 +- configure.ac | 2 +- ltmain.sh | 4 +- m4/libtool.m4 | 12 +- src/algorithms/ImgRegression.cc | 309 +++++++++++++- src/algorithms/ImgRegression.h | 8 +- src/algorithms/StatFactory.h | 33 +- src/apps/Makefile.am | 6 +- src/apps/Makefile.in | 71 ++-- src/apps/pkann.cc | 72 ++-- src/apps/pkcomposite.cc | 62 +-- src/apps/pkcrop.cc | 7 +- src/apps/pkdiff.cc | 29 +- src/apps/pkdumpimg.cc | 14 +- src/apps/pkenhance.cc | 4 +- src/apps/pkextract.cc | 4 - src/apps/pkinfo.cc | 53 +-- src/apps/pkkalman.cc | 603 ++++++++++++++++++++-------- src/apps/pkstat.cc | 825 ++++++++++++++++++++++++++++++++++++++ src/apps/pkstatascii.cc | 12 +- src/apps/pkstatogr.cc | 2 +- src/apps/pksvm.cc | 70 ++-- src/imageclasses/ImgReaderGdal.cc | 73 +++- src/imageclasses/ImgReaderGdal.h | 3 +- src/imageclasses/Makefile.am | 2 +- src/imageclasses/Makefile.in | 2 +- 27 files changed, 1904 insertions(+), 426 deletions(-) diff --git a/ChangeLog b/ChangeLog index 59a04c0..dad167d 100755 --- a/ChangeLog +++ b/ChangeLog @@ -346,6 +346,22 @@ version 2.6.2 support of duplicate option -b to select band in reference image always use first band in mask image instead of using band defined in option -b - command line help info (change request #43845) + +version 2.6.3 + - pksvm + solved bug in mask handling + - pkann + solved bug in mask handling + - pkcomposite + Advanced option -file supports two modes: 1 (total number of files used) and 2 (sequence number of selected file) + No default for option -srcnodata + - pkstat + New utility for calculating statistics on raster datasets (regression, histograms, etc.) + - pkinfo + Removed hist from pkinfo (now in pkstat) + - pkdiff + support double data types for input, but cast to unsigned short in case of accuracy assessment + Todo: - todo for API ImgReaderGdal (ImgWriterGdal) open in update mode (check gdal_edit.py: http://searchcode.com/codesearch/view/18938404) diff --git a/configure b/configure index 9106467..3ab8c93 100755 --- a/configure +++ b/configure @@ -1,6 +1,6 @@ #! /bin/sh # Guess values for system-dependent variables and create Makefiles. -# Generated by GNU Autoconf 2.69 for pktools 2.6.2. +# Generated by GNU Autoconf 2.69 for pktools 2.6.3. # # Report bugs to <[email protected]>. # @@ -590,8 +590,8 @@ MAKEFLAGS= # Identity of this package. PACKAGE_NAME='pktools' PACKAGE_TARNAME='pktools' -PACKAGE_VERSION='2.6.2' -PACKAGE_STRING='pktools 2.6.2' +PACKAGE_VERSION='2.6.3' +PACKAGE_STRING='pktools 2.6.3' PACKAGE_BUGREPORT='[email protected]' PACKAGE_URL='' @@ -1366,7 +1366,7 @@ if test "$ac_init_help" = "long"; then # Omit some internal or obsolete options to make the list less imposing. # This message is too long to be a string in the A/UX 3.1 sh. cat <<_ACEOF -\`configure' configures pktools 2.6.2 to adapt to many kinds of systems. +\`configure' configures pktools 2.6.3 to adapt to many kinds of systems. Usage: $0 [OPTION]... [VAR=VALUE]... @@ -1436,7 +1436,7 @@ fi if test -n "$ac_init_help"; then case $ac_init_help in - short | recursive ) echo "Configuration of pktools 2.6.2:";; + short | recursive ) echo "Configuration of pktools 2.6.3:";; esac cat <<\_ACEOF @@ -1562,7 +1562,7 @@ fi test -n "$ac_init_help" && exit $ac_status if $ac_init_version; then cat <<\_ACEOF -pktools configure 2.6.2 +pktools configure 2.6.3 generated by GNU Autoconf 2.69 Copyright (C) 2012 Free Software Foundation, Inc. @@ -2323,7 +2323,7 @@ cat >config.log <<_ACEOF This file contains any messages produced by compilers while running configure, to aid debugging if configure makes a mistake. -It was created by pktools $as_me 2.6.2, which was +It was created by pktools $as_me 2.6.3, which was generated by GNU Autoconf 2.69. Invocation command line was $ $0 $@ @@ -3187,7 +3187,7 @@ fi # Define the identity of the package. PACKAGE='pktools' - VERSION='2.6.2' + VERSION='2.6.3' cat >>confdefs.h <<_ACEOF @@ -7520,7 +7520,7 @@ ia64-*-hpux*) rm -rf conftest* ;; -x86_64-*kfreebsd*-gnu|x86_64-*linux*|powerpc*-*linux*| \ +x86_64-*kfreebsd*-gnu|x86_64-*linux*|ppc*-*linux*|powerpc*-*linux*| \ s390*-*linux*|s390*-*tpf*|sparc*-*linux*) # Find out which ABI we are using. echo 'int i;' > conftest.$ac_ext @@ -7545,10 +7545,7 @@ s390*-*linux*|s390*-*tpf*|sparc*-*linux*) ;; esac ;; - powerpc64le-*) - LD="${LD-ld} -m elf32lppclinux" - ;; - powerpc64-*) + ppc64-*linux*|powerpc64-*linux*) LD="${LD-ld} -m elf32ppclinux" ;; s390x-*linux*) @@ -7567,10 +7564,7 @@ s390*-*linux*|s390*-*tpf*|sparc*-*linux*) x86_64-*linux*) LD="${LD-ld} -m elf_x86_64" ;; - powerpcle-*) - LD="${LD-ld} -m elf64lppc" - ;; - powerpc-*) + ppc*-*linux*|powerpc*-*linux*) LD="${LD-ld} -m elf64ppc" ;; s390*-*linux*|s390*-*tpf*) @@ -20172,7 +20166,7 @@ cat >>$CONFIG_STATUS <<\_ACEOF || ac_write_fail=1 # report actual input values of CONFIG_FILES etc. instead of their # values after options handling. ac_log=" -This file was extended by pktools $as_me 2.6.2, which was +This file was extended by pktools $as_me 2.6.3, which was generated by GNU Autoconf 2.69. Invocation command line was CONFIG_FILES = $CONFIG_FILES @@ -20238,7 +20232,7 @@ _ACEOF cat >>$CONFIG_STATUS <<_ACEOF || ac_write_fail=1 ac_cs_config="`$as_echo "$ac_configure_args" | sed 's/^ //; s/[\\""\`\$]/\\\\&/g'`" ac_cs_version="\\ -pktools config.status 2.6.2 +pktools config.status 2.6.3 configured by $0, generated by GNU Autoconf 2.69, with options \\"\$ac_cs_config\\" diff --git a/configure.ac b/configure.ac index 631ff56..98bef35 100644 --- a/configure.ac +++ b/configure.ac @@ -1,4 +1,4 @@ -AC_INIT([pktools], [2.6.2], [[email protected]]) +AC_INIT([pktools], [2.6.3], [[email protected]]) #AM_INIT_AUTOMAKE([-Wall -Werror foreign]) AM_INIT_AUTOMAKE([-Wall -Wno-extra-portability foreign]) #AM_INIT_AUTOMAKE([subdir-objects]) #not working due to bug in autoconf, see Debian list: Bug #752993) diff --git a/ltmain.sh b/ltmain.sh index a356aca..b9205ee 100644 --- a/ltmain.sh +++ b/ltmain.sh @@ -70,7 +70,7 @@ # compiler: $LTCC # compiler flags: $LTCFLAGS # linker: $LD (gnu? $with_gnu_ld) -# $progname: (GNU libtool) 2.4.2 Debian-2.4.2-1.7ubuntu1 +# $progname: (GNU libtool) 2.4.2 Debian-2.4.2-1.2ubuntu1 # automake: $automake_version # autoconf: $autoconf_version # @@ -80,7 +80,7 @@ PROGRAM=libtool PACKAGE=libtool -VERSION="2.4.2 Debian-2.4.2-1.7ubuntu1" +VERSION="2.4.2 Debian-2.4.2-1.2ubuntu1" TIMESTAMP="" package_revision=1.3337 diff --git a/m4/libtool.m4 b/m4/libtool.m4 index d7c043f..02b4bbe 100644 --- a/m4/libtool.m4 +++ b/m4/libtool.m4 @@ -1312,7 +1312,7 @@ ia64-*-hpux*) rm -rf conftest* ;; -x86_64-*kfreebsd*-gnu|x86_64-*linux*|powerpc*-*linux*| \ +x86_64-*kfreebsd*-gnu|x86_64-*linux*|ppc*-*linux*|powerpc*-*linux*| \ s390*-*linux*|s390*-*tpf*|sparc*-*linux*) # Find out which ABI we are using. echo 'int i;' > conftest.$ac_ext @@ -1333,10 +1333,7 @@ s390*-*linux*|s390*-*tpf*|sparc*-*linux*) ;; esac ;; - powerpc64le-*) - LD="${LD-ld} -m elf32lppclinux" - ;; - powerpc64-*) + ppc64-*linux*|powerpc64-*linux*) LD="${LD-ld} -m elf32ppclinux" ;; s390x-*linux*) @@ -1355,10 +1352,7 @@ s390*-*linux*|s390*-*tpf*|sparc*-*linux*) x86_64-*linux*) LD="${LD-ld} -m elf_x86_64" ;; - powerpcle-*) - LD="${LD-ld} -m elf64lppc" - ;; - powerpc-*) + ppc*-*linux*|powerpc*-*linux*) LD="${LD-ld} -m elf64ppc" ;; s390*-*linux*|s390*-*tpf*) diff --git a/src/algorithms/ImgRegression.cc b/src/algorithms/ImgRegression.cc index e1de847..6b4d727 100644 --- a/src/algorithms/ImgRegression.cc +++ b/src/algorithms/ImgRegression.cc @@ -29,7 +29,7 @@ ImgRegression::ImgRegression(void) ImgRegression::~ImgRegression(void) {} -double ImgRegression::getRMSE(const ImgReaderGdal& imgReader1, const ImgReaderGdal& imgReader2, double& c0, double& c1, short verbose) const{ +double ImgRegression::getRMSE(const ImgReaderGdal& imgReader1, const ImgReaderGdal& imgReader2, double& c0, double& c1, unsigned short band1, unsigned short band2, short verbose) const{ c0=0; c1=1; int icol1=0,irow1=0; @@ -45,14 +45,14 @@ double ImgRegression::getRMSE(const ImgReaderGdal& imgReader1, const ImgReaderGd icol1=0; double icol2=0,irow2=0; double geox=0,geoy=0; - imgReader1.readData(rowBuffer1,GDT_Float64,irow1); + imgReader1.readData(rowBuffer1,GDT_Float64,irow1,band1); imgReader1.image2geo(icol1,irow1,geox,geoy); imgReader2.geo2image(geox,geoy,icol2,irow2); icol2=static_cast<int>(icol2); irow2=static_cast<int>(irow2); if(irow2<0||irow2>=imgReader2.nrOfRow()) continue; - imgReader2.readData(rowBuffer2,GDT_Float64,irow2); + imgReader2.readData(rowBuffer2,GDT_Float64,irow2,band2); for(icol1=0;icol1<imgReader1.nrOfCol();++icol1){ if(icol1%m_down) continue; @@ -89,3 +89,306 @@ double ImgRegression::getRMSE(const ImgReaderGdal& imgReader1, const ImgReaderGd std::cout << "linear regression based on " << buffer1.size() << " points: " << c0 << "+" << c1 << " * x " << " with rmse: " << err << std::endl; return err; } + +double ImgRegression::getR2(const ImgReaderGdal& imgReader1, const ImgReaderGdal& imgReader2, double& c0, double& c1, unsigned short band1, unsigned short band2, short verbose) const{ + c0=0; + c1=1; + int icol1=0,irow1=0; + std::vector<double> rowBuffer1(imgReader1.nrOfCol()); + std::vector<double> rowBuffer2(imgReader2.nrOfCol()); + std::vector<double> buffer1; + std::vector<double> buffer2; + + srand(time(NULL)); + for(irow1=0;irow1<imgReader1.nrOfRow();++irow1){ + if(irow1%m_down) + continue; + icol1=0; + double icol2=0,irow2=0; + double geox=0,geoy=0; + imgReader1.readData(rowBuffer1,GDT_Float64,irow1,band1); + imgReader1.image2geo(icol1,irow1,geox,geoy); + imgReader2.geo2image(geox,geoy,icol2,irow2); + icol2=static_cast<int>(icol2); + irow2=static_cast<int>(irow2); + if(irow2<0||irow2>=imgReader2.nrOfRow()) + continue; + imgReader2.readData(rowBuffer2,GDT_Float64,irow2,band2); + for(icol1=0;icol1<imgReader1.nrOfCol();++icol1){ + if(icol1%m_down) + continue; + if(m_threshold>0){//percentual value + double p=static_cast<double>(rand())/(RAND_MAX); + p*=100.0; + if(p>m_threshold) + continue;//do not select for now, go to next column + } + imgReader1.image2geo(icol1,irow1,geox,geoy); + imgReader2.geo2image(geox,geoy,icol2,irow2); + if(icol2<0||icol2>=imgReader2.nrOfCol()) + continue; + icol2=static_cast<int>(icol2); + irow2=static_cast<int>(irow2); + //check for nodata + double value1=rowBuffer1[icol1]; + double value2=rowBuffer2[icol2]; + if(imgReader1.isNoData(value1)||imgReader2.isNoData(value2)) + continue; + + buffer1.push_back(value1); + buffer2.push_back(value2); + if(verbose>1) + std::cout << geox << " " << geoy << " " << icol1 << " " << irow1 << " " << icol2 << " " << irow2 << " " << buffer1.back() << " " << buffer2.back() << std::endl; + } + } + double r2=0; + if(buffer1.size()&&buffer2.size()){ + statfactory::StatFactory stat; + r2=stat.linear_regression(buffer1,buffer2,c0,c1); + } + if(verbose) + std::cout << "linear regression based on " << buffer1.size() << " points: " << c0 << "+" << c1 << " * x " << " with r^2: " << r2 << std::endl; + return r2; +} + +double ImgRegression::pgetR2(const ImgReaderGdal& imgReader1, const ImgReaderGdal& imgReader2, double& c0, double& c1, unsigned short band1, unsigned short band2, short verbose) const{ + c0=0; + c1=1; + int icol1=0,irow1=0; + std::vector<double> rowBuffer1(imgReader1.nrOfCol()); + std::vector<double> rowBuffer2(imgReader2.nrOfCol()); + std::vector<double> buffer1; + std::vector<double> buffer2; + + srand(time(NULL)); + for(irow1=0;irow1<imgReader1.nrOfRow();++irow1){ + if(irow1%m_down) + continue; + icol1=0; + double icol2=0,irow2=0; + double geox=0,geoy=0; + imgReader1.readData(rowBuffer1,GDT_Float64,irow1,band1); + imgReader1.image2geo(icol1,irow1,geox,geoy); + imgReader2.geo2image(geox,geoy,icol2,irow2); + icol2=static_cast<int>(icol2); + irow2=static_cast<int>(irow2); + if(irow2<0||irow2>=imgReader2.nrOfRow()) + continue; + imgReader2.readData(rowBuffer2,GDT_Float64,irow2,band2); + for(icol1=0;icol1<imgReader1.nrOfCol();++icol1){ + if(icol1%m_down) + continue; + if(m_threshold>0){//percentual value + double p=static_cast<double>(rand())/(RAND_MAX); + p*=100.0; + if(p>m_threshold) + continue;//do not select for now, go to next column + } + imgReader1.image2geo(icol1,irow1,geox,geoy); + imgReader2.geo2image(geox,geoy,icol2,irow2); + if(icol2<0||icol2>=imgReader2.nrOfCol()) + continue; + icol2=static_cast<int>(icol2); + irow2=static_cast<int>(irow2); + //check for nodata + double value1=rowBuffer1[icol1]; + double value2=rowBuffer2[icol2]; + if(imgReader1.isNoData(value1)||imgReader2.isNoData(value2)) + continue; + + buffer1.push_back(value1); + buffer2.push_back(value2); + if(verbose>1) + std::cout << geox << " " << geoy << " " << icol1 << " " << irow1 << " " << icol2 << " " << irow2 << " " << buffer1.back() << " " << buffer2.back() << std::endl; + } + } + double r=0; + if(buffer1.size()&&buffer2.size()){ + statfactory::StatFactory stat; + r=stat.correlation(buffer1,buffer2); + // r=stat.gsl_correlation(buffer1,buffer2); + double m1=0; + double v1=0; + double m2=0; + double v2=0; + stat.meanVar(buffer1,m1,v1); + stat.meanVar(buffer2,m2,v2); + if(v1>0){ + if(r>=0) + c1=v2/v1; + else + c1=-v2/v1; + } + c0=m2-c1*m1; + } + if(verbose) + std::cout << "orthogonal regression based on " << buffer1.size() << " points: " << c0 << "+" << c1 << " * x " << " with r^2: " << r*r << std::endl; + return r*r; +} + +double ImgRegression::getRMSE(const ImgReaderGdal& imgReader, unsigned short band1, unsigned short band2, double& c0, double& c1, short verbose) const{ + c0=0; + c1=1; + int icol=0,irow=0; + std::vector<double> rowBuffer1(imgReader.nrOfCol()); + std::vector<double> rowBuffer2(imgReader.nrOfCol()); + std::vector<double> buffer1; + std::vector<double> buffer2; + + srand(time(NULL)); + assert(band1>=0); + assert(band1<imgReader.nrOfBand()); + assert(band2>=0); + assert(band2<imgReader.nrOfBand()); + for(irow=0;irow<imgReader.nrOfRow();++irow){ + if(irow%m_down) + continue; + icol=0; + imgReader.readData(rowBuffer1,GDT_Float64,irow,band1); + imgReader.readData(rowBuffer2,GDT_Float64,irow,band2); + for(icol=0;icol<imgReader.nrOfCol();++icol){ + if(icol%m_down) + continue; + if(m_threshold>0){//percentual value + double p=static_cast<double>(rand())/(RAND_MAX); + p*=100.0; + if(p>m_threshold) + continue;//do not select for now, go to next column + } + //check for nodata + double value1=rowBuffer1[icol]; + double value2=rowBuffer2[icol]; + if(imgReader.isNoData(value1)||imgReader.isNoData(value2)) + continue; + + buffer1.push_back(value1); + buffer2.push_back(value2); + if(verbose>1) + std::cout << icol << " " << irow << " " << buffer1.back() << " " << buffer2.back() << std::endl; + } + } + double err=0; + if(buffer1.size()&&buffer2.size()){ + statfactory::StatFactory stat; + err=stat.linear_regression_err(buffer1,buffer2,c0,c1); + } + if(verbose) + std::cout << "linear regression based on " << buffer1.size() << " points: " << c0 << "+" << c1 << " * x " << " with rmse: " << err << std::endl; + return err; +} + +double ImgRegression::getR2(const ImgReaderGdal& imgReader, unsigned short band1, unsigned short band2, double& c0, double& c1, short verbose) const{ + c0=0; + c1=1; + int icol=0,irow=0; + std::vector<double> rowBuffer1(imgReader.nrOfCol()); + std::vector<double> rowBuffer2(imgReader.nrOfCol()); + std::vector<double> buffer1; + std::vector<double> buffer2; + + srand(time(NULL)); + assert(band1>=0); + assert(band1<imgReader.nrOfBand()); + assert(band2>=0); + assert(band2<imgReader.nrOfBand()); + for(irow=0;irow<imgReader.nrOfRow();++irow){ + if(irow%m_down) + continue; + icol=0; + imgReader.readData(rowBuffer1,GDT_Float64,irow,band1); + imgReader.readData(rowBuffer2,GDT_Float64,irow,band2); + for(icol=0;icol<imgReader.nrOfCol();++icol){ + if(icol%m_down) + continue; + if(m_threshold>0){//percentual value + double p=static_cast<double>(rand())/(RAND_MAX); + p*=100.0; + if(p>m_threshold) + continue;//do not select for now, go to next column + } + //check for nodata + double value1=rowBuffer1[icol]; + double value2=rowBuffer2[icol]; + if(imgReader.isNoData(value1)||imgReader.isNoData(value2)) + continue; + + buffer1.push_back(value1); + buffer2.push_back(value2); + if(verbose>1) + std::cout << icol << " " << irow << " " << buffer1.back() << " " << buffer2.back() << std::endl; + } + } + double r2=0; + if(buffer1.size()&&buffer2.size()){ + statfactory::StatFactory stat; + r2=stat.linear_regression(buffer1,buffer2,c0,c1); + } + if(verbose) + std::cout << "linear regression based on " << buffer1.size() << " points: " << c0 << "+" << c1 << " * x " << " with r^2: " << r2 << std::endl; + return r2; +} + +double ImgRegression::pgetR2(const ImgReaderGdal& imgReader, unsigned short band1, unsigned short band2, double& c0, double& c1, short verbose) const{ + c0=0; + c1=1; + int icol=0,irow=0; + std::vector<double> rowBuffer1(imgReader.nrOfCol()); + std::vector<double> rowBuffer2(imgReader.nrOfCol()); + std::vector<double> buffer1; + std::vector<double> buffer2; + + srand(time(NULL)); + assert(band1>=0); + assert(band1<imgReader.nrOfBand()); + assert(band2>=0); + assert(band2<imgReader.nrOfBand()); + for(irow=0;irow<imgReader.nrOfRow();++irow){ + if(irow%m_down) + continue; + icol=0; + imgReader.readData(rowBuffer1,GDT_Float64,irow,band1); + imgReader.readData(rowBuffer2,GDT_Float64,irow,band2); + for(icol=0;icol<imgReader.nrOfCol();++icol){ + if(icol%m_down) + continue; + if(m_threshold>0){//percentual value + double p=static_cast<double>(rand())/(RAND_MAX); + p*=100.0; + if(p>m_threshold) + continue;//do not select for now, go to next column + } + //check for nodata + double value1=rowBuffer1[icol]; + double value2=rowBuffer2[icol]; + if(imgReader.isNoData(value1)||imgReader.isNoData(value2)) + continue; + + buffer1.push_back(value1); + buffer2.push_back(value2); + if(verbose>1) + std::cout << icol << " " << irow << " " << buffer1.back() << " " << buffer2.back() << std::endl; + } + } + double r=0; + if(buffer1.size()&&buffer2.size()){ + statfactory::StatFactory stat; + r=stat.correlation(buffer1,buffer2); + // r=stat.gsl_correlation(buffer1,buffer2); + double m1=0; + double v1=0; + double m2=0; + double v2=0; + stat.meanVar(buffer1,m1,v1); + stat.meanVar(buffer2,m2,v2); + if(v1>0){ + if(r>=0) + c1=v2/v1; + else + c1=-v2/v1; + } + c0=m2-c1*m1; + } + if(verbose) + std::cout << "orthogonal regression based on " << buffer1.size() << " points: " << c0 << "+" << c1 << " * x " << " with r^2: " << r*r << std::endl; + return r*r; +} diff --git a/src/algorithms/ImgRegression.h b/src/algorithms/ImgRegression.h index 12ba62d..51bea30 100644 --- a/src/algorithms/ImgRegression.h +++ b/src/algorithms/ImgRegression.h @@ -31,7 +31,13 @@ namespace imgregression public: ImgRegression(void); ~ImgRegression(void); - double getRMSE(const ImgReaderGdal& imgReader1, const ImgReaderGdal& imgReader2, double &c0, double &c1, short verbose=0) const; + double getRMSE(const ImgReaderGdal& imgReader1, const ImgReaderGdal& imgReader2, double &c0, double &c1, unsigned short b1=0, unsigned short b2=0, short verbose=0) const; + double getRMSE(const ImgReaderGdal& imgReader, unsigned short b1, unsigned short b2, double& c0, double& c1, short verbose=0) const; + double getR2(const ImgReaderGdal& imgReader1, const ImgReaderGdal& imgReader2, double &c0, double &c1, unsigned short b1=0, unsigned short b2=0, short verbose=0) const; + double pgetR2(const ImgReaderGdal& imgReader1, const ImgReaderGdal& imgReader2, double& c0, double& c1, unsigned short band1, unsigned short band2, short verbose=0) const; + double getR2(const ImgReaderGdal& imgReader, unsigned short b1, unsigned short b2, double& c0, double& c1, short verbose=0) const; + double pgetR2(const ImgReaderGdal& imgReader, unsigned short band1, unsigned short band2, double& c0, double& c1, short verbose=0) const; + void setThreshold(double theThreshold){m_threshold=theThreshold;}; void setDown(int theDown){m_down=theDown;}; private: diff --git a/src/algorithms/StatFactory.h b/src/algorithms/StatFactory.h index b9511d8..b307486 100644 --- a/src/algorithms/StatFactory.h +++ b/src/algorithms/StatFactory.h @@ -176,6 +176,8 @@ public: template<class T> void normalize_pct(std::vector<T>& input) const; template<class T> double rmse(const std::vector<T>& x, const std::vector<T>& y) const; template<class T> double correlation(const std::vector<T>& x, const std::vector<T>& y, int delay=0) const; + // template<class T> double gsl_correlation(const std::vector<T>& x, const std::vector<T>& y) const; + template<class T> double gsl_covariance(const std::vector<T>& x, const std::vector<T>& y) const; template<class T> double cross_correlation(const std::vector<T>& x, const std::vector<T>& y, int maxdelay, std::vector<T>& z) const; template<class T> double linear_regression(const std::vector<T>& x, const std::vector<T>& y, double &c0, double &c1) const; template<class T> double linear_regression_err(const std::vector<T>& x, const std::vector<T>& y, double &c0, double &c1) const; @@ -875,10 +877,24 @@ template<class T> void StatFactory::distribution2d(const std::vector<T>& inputX s<<"Error opening distribution file , " << filename; throw(s.str()); } - for(int bin=0;bin<nbin;++bin){ - for(int bin=0;bin<nbin;++bin){ - double value=static_cast<double>(output[binX][binY])/npoint; - outputfile << (maxX-minX)*bin/(nbin-1)+minX << " " << (maxY-minY)*bin/(nbin-1)+minY << " " << value << std::endl; + for(int binX=0;binX<nbin;++binX){ + outputfile << std::endl; + for(int binY=0;binY<nbin;++binY){ + double binValueX=0; + if(nbin==maxX-minX+1) + binValueX=minX+binX; + else + binValueX=minX+static_cast<double>(maxX-minX)*(binX+0.5)/nbin; + double binValueY=0; + if(nbin==maxY-minY+1) + binValueY=minY+binY; + else + binValueY=minY+static_cast<double>(maxY-minY)*(binY+0.5)/nbin; + double value=0; + value=static_cast<double>(output[binX][binY])/npoint; + outputfile << binValueX << " " << binValueY << " " << value << std::endl; + /* double value=static_cast<double>(output[binX][binY])/npoint; */ + /* outputfile << (maxX-minX)*bin/(nbin-1)+minX << " " << (maxY-minY)*bin/(nbin-1)+minY << " " << value << std::endl; */ } } outputfile.close(); @@ -972,6 +988,15 @@ template<class T> double StatFactory::rmse(const std::vector<T>& x, const std::v return sqrt(mse); } +// template<class T> double StatFactory::gsl_correlation(const std::vector<T>& x, const std::vector<T>& y) const{ +// return(gsl_stats_correlation(&(x[0]),1,&(y[0]),1,x.size())); +// } + + template<class T> double StatFactory::gsl_covariance(const std::vector<T>& x, const std::vector<T>& y) const{ + return(gsl_stats_covariance(&(x[0]),1,&(y[0]),1,x.size())); + } + + template<class T> double StatFactory::correlation(const std::vector<T>& x, const std::vector<T>& y, int delay) const{ double meanX=0; double meanY=0; diff --git a/src/apps/Makefile.am b/src/apps/Makefile.am index a2a243a..c3f3be3 100644 --- a/src/apps/Makefile.am +++ b/src/apps/Makefile.am @@ -6,10 +6,10 @@ LDADD = $(GSL_LIBS) $(GDAL_LDFLAGS) $(top_builddir)/src/algorithms/libalgorithms ############################################################################### # the program to build and install (the names of the final binaries) -bin_PROGRAMS = pkinfo pkcrop pkreclass pkdiff pkgetmask pksetmask pkcreatect pkdumpimg pkdumpogr pksieve pkstatascii pkstatogr pkegcs pkextract pkfillnodata pkfilter pkkalman pkfilterdem pkenhance pkfilterascii pkdsm2shadow pkcomposite pkndvi pkpolygonize pkascii2img pksvm pkfssvm pkascii2ogr pkeditogr +bin_PROGRAMS = pkinfo pkcrop pkdiff pkgetmask pksetmask pkcreatect pkdumpimg pkdumpogr pksieve pkstat pkstatascii pkstatogr pkegcs pkextract pkfillnodata pkfilter pkfilterdem pkfilterascii pkdsm2shadow pkcomposite pkpolygonize pkascii2img pksvm pkfssvm pkascii2ogr pkkalman # the program to build but not install (the names of the final binaries) -#noinst_PROGRAMS = pkxcorimg pkgeom +noinst_PROGRAMS = pkeditogr pkenhance pkndvi pkreclass if USE_FANN bin_PROGRAMS += pkann pkfsann pkregann @@ -48,6 +48,8 @@ pkcreatect_SOURCES = pkcreatect.cc pkdumpimg_SOURCES = pkdumpimg.cc pkdumpogr_SOURCES = pkdumpogr.h pkdumpogr.cc pksieve_SOURCES = pksieve.cc +pkstat_SOURCES = $(top_srcdir)/src/algorithms/StatFactory.h $(top_srcdir)/src/algorithms/ImgRegression.h $(top_srcdir)/src/algorithms/ImgRegression.cc pkstat.cc +pkstat_LDADD = -lgsl -lgslcblas $(GSL_LIBS) $(AM_LDFLAGS) pkstatascii_SOURCES = $(top_srcdir)/src/algorithms/StatFactory.h pkstatascii.cc pkstatascii_LDADD = -lgsl -lgslcblas $(GSL_LIBS) $(AM_LDFLAGS) pkstatogr_SOURCES = pkstatogr.cc diff --git a/src/apps/Makefile.in b/src/apps/Makefile.in index 363a78f..90adb9b 100644 --- a/src/apps/Makefile.in +++ b/src/apps/Makefile.in @@ -78,20 +78,18 @@ PRE_UNINSTALL = : POST_UNINSTALL = : build_triplet = @build@ host_triplet = @host@ -bin_PROGRAMS = pkinfo$(EXEEXT) pkcrop$(EXEEXT) pkreclass$(EXEEXT) \ - pkdiff$(EXEEXT) pkgetmask$(EXEEXT) pksetmask$(EXEEXT) \ - pkcreatect$(EXEEXT) pkdumpimg$(EXEEXT) pkdumpogr$(EXEEXT) \ - pksieve$(EXEEXT) pkstatascii$(EXEEXT) pkstatogr$(EXEEXT) \ +bin_PROGRAMS = pkinfo$(EXEEXT) pkcrop$(EXEEXT) pkdiff$(EXEEXT) \ + pkgetmask$(EXEEXT) pksetmask$(EXEEXT) pkcreatect$(EXEEXT) \ + pkdumpimg$(EXEEXT) pkdumpogr$(EXEEXT) pksieve$(EXEEXT) \ + pkstat$(EXEEXT) pkstatascii$(EXEEXT) pkstatogr$(EXEEXT) \ pkegcs$(EXEEXT) pkextract$(EXEEXT) pkfillnodata$(EXEEXT) \ - pkfilter$(EXEEXT) pkkalman$(EXEEXT) pkfilterdem$(EXEEXT) \ - pkenhance$(EXEEXT) pkfilterascii$(EXEEXT) \ - pkdsm2shadow$(EXEEXT) pkcomposite$(EXEEXT) pkndvi$(EXEEXT) \ + pkfilter$(EXEEXT) pkfilterdem$(EXEEXT) pkfilterascii$(EXEEXT) \ + pkdsm2shadow$(EXEEXT) pkcomposite$(EXEEXT) \ pkpolygonize$(EXEEXT) pkascii2img$(EXEEXT) pksvm$(EXEEXT) \ - pkfssvm$(EXEEXT) pkascii2ogr$(EXEEXT) pkeditogr$(EXEEXT) \ + pkfssvm$(EXEEXT) pkascii2ogr$(EXEEXT) pkkalman$(EXEEXT) \ $(am__EXEEXT_1) $(am__EXEEXT_2) $(am__EXEEXT_3) - -# the program to build but not install (the names of the final binaries) -#noinst_PROGRAMS = pkxcorimg pkgeom +noinst_PROGRAMS = pkeditogr$(EXEEXT) pkenhance$(EXEEXT) \ + pkndvi$(EXEEXT) pkreclass$(EXEEXT) @USE_FANN_TRUE@am__append_1 = pkann pkfsann pkregann @USE_LAS_TRUE@am__append_2 = pklas2img @USE_NLOPT_TRUE@am__append_3 = pkoptsvm @@ -114,7 +112,7 @@ CONFIG_CLEAN_VPATH_FILES = @USE_LAS_TRUE@am__EXEEXT_2 = pklas2img$(EXEEXT) @USE_NLOPT_TRUE@am__EXEEXT_3 = pkoptsvm$(EXEEXT) am__installdirs = "$(DESTDIR)$(bindir)" -PROGRAMS = $(bin_PROGRAMS) +PROGRAMS = $(bin_PROGRAMS) $(noinst_PROGRAMS) am__pkann_SOURCES_DIST = $(top_srcdir)/src/algorithms/myfann_cpp.h \ pkann.cc @USE_FANN_TRUE@am_pkann_OBJECTS = pkann-pkann.$(OBJEXT) @@ -303,6 +301,9 @@ pksieve_DEPENDENCIES = $(am__DEPENDENCIES_1) $(am__DEPENDENCIES_1) \ $(top_builddir)/src/algorithms/libalgorithms.la \ $(top_builddir)/src/imageclasses/libimageClasses.la \ $(top_builddir)/src/fileclasses/libfileClasses.la +am_pkstat_OBJECTS = ImgRegression.$(OBJEXT) pkstat.$(OBJEXT) +pkstat_OBJECTS = $(am_pkstat_OBJECTS) +pkstat_DEPENDENCIES = $(am__DEPENDENCIES_1) $(am__DEPENDENCIES_2) am_pkstatascii_OBJECTS = pkstatascii.$(OBJEXT) pkstatascii_OBJECTS = $(am_pkstatascii_OBJECTS) pkstatascii_DEPENDENCIES = $(am__DEPENDENCIES_1) $(am__DEPENDENCIES_2) @@ -376,8 +377,8 @@ SOURCES = $(pkann_SOURCES) $(pkascii2img_SOURCES) \ $(pkinfo_SOURCES) $(pkkalman_SOURCES) $(pklas2img_SOURCES) \ $(pkndvi_SOURCES) $(pkoptsvm_SOURCES) $(pkpolygonize_SOURCES) \ $(pkreclass_SOURCES) $(pkregann_SOURCES) $(pksetmask_SOURCES) \ - $(pksieve_SOURCES) $(pkstatascii_SOURCES) $(pkstatogr_SOURCES) \ - $(pksvm_SOURCES) + $(pksieve_SOURCES) $(pkstat_SOURCES) $(pkstatascii_SOURCES) \ + $(pkstatogr_SOURCES) $(pksvm_SOURCES) DIST_SOURCES = $(am__pkann_SOURCES_DIST) $(pkascii2img_SOURCES) \ $(pkascii2ogr_SOURCES) $(pkcomposite_SOURCES) \ $(pkcreatect_SOURCES) $(pkcrop_SOURCES) $(pkdiff_SOURCES) \ @@ -391,8 +392,8 @@ DIST_SOURCES = $(am__pkann_SOURCES_DIST) $(pkascii2img_SOURCES) \ $(am__pklas2img_SOURCES_DIST) $(pkndvi_SOURCES) \ $(am__pkoptsvm_SOURCES_DIST) $(pkpolygonize_SOURCES) \ $(pkreclass_SOURCES) $(am__pkregann_SOURCES_DIST) \ - $(pksetmask_SOURCES) $(pksieve_SOURCES) $(pkstatascii_SOURCES) \ - $(pkstatogr_SOURCES) $(pksvm_SOURCES) + $(pksetmask_SOURCES) $(pksieve_SOURCES) $(pkstat_SOURCES) \ + $(pkstatascii_SOURCES) $(pkstatogr_SOURCES) $(pksvm_SOURCES) am__can_run_installinfo = \ case $$AM_UPDATE_INFO_DIR in \ n|no|NO) false;; \ @@ -587,6 +588,8 @@ pkcreatect_SOURCES = pkcreatect.cc pkdumpimg_SOURCES = pkdumpimg.cc pkdumpogr_SOURCES = pkdumpogr.h pkdumpogr.cc pksieve_SOURCES = pksieve.cc +pkstat_SOURCES = $(top_srcdir)/src/algorithms/StatFactory.h $(top_srcdir)/src/algorithms/ImgRegression.h $(top_srcdir)/src/algorithms/ImgRegression.cc pkstat.cc +pkstat_LDADD = -lgsl -lgslcblas $(GSL_LIBS) $(AM_LDFLAGS) pkstatascii_SOURCES = $(top_srcdir)/src/algorithms/StatFactory.h pkstatascii.cc pkstatascii_LDADD = -lgsl -lgslcblas $(GSL_LIBS) $(AM_LDFLAGS) pkstatogr_SOURCES = pkstatogr.cc @@ -702,6 +705,15 @@ clean-binPROGRAMS: echo " rm -f" $$list; \ rm -f $$list +clean-noinstPROGRAMS: + @list='$(noinst_PROGRAMS)'; test -n "$$list" || exit 0; \ + echo " rm -f" $$list; \ + rm -f $$list || exit $$?; \ + test -n "$(EXEEXT)" || exit 0; \ + list=`for p in $$list; do echo "$$p"; done | sed 's/$(EXEEXT)$$//'`; \ + echo " rm -f" $$list; \ + rm -f $$list + pkann$(EXEEXT): $(pkann_OBJECTS) $(pkann_DEPENDENCIES) $(EXTRA_pkann_DEPENDENCIES) @rm -f pkann$(EXEEXT) $(AM_V_CXXLD)$(pkann_LINK) $(pkann_OBJECTS) $(pkann_LDADD) $(LIBS) @@ -826,6 +838,10 @@ pksieve$(EXEEXT): $(pksieve_OBJECTS) $(pksieve_DEPENDENCIES) $(EXTRA_pksieve_DEP @rm -f pksieve$(EXEEXT) $(AM_V_CXXLD)$(CXXLINK) $(pksieve_OBJECTS) $(pksieve_LDADD) $(LIBS) +pkstat$(EXEEXT): $(pkstat_OBJECTS) $(pkstat_DEPENDENCIES) $(EXTRA_pkstat_DEPENDENCIES) + @rm -f pkstat$(EXEEXT) + $(AM_V_CXXLD)$(CXXLINK) $(pkstat_OBJECTS) $(pkstat_LDADD) $(LIBS) + pkstatascii$(EXEEXT): $(pkstatascii_OBJECTS) $(pkstatascii_DEPENDENCIES) $(EXTRA_pkstatascii_DEPENDENCIES) @rm -f pkstatascii$(EXEEXT) $(AM_V_CXXLD)$(CXXLINK) $(pkstatascii_OBJECTS) $(pkstatascii_LDADD) $(LIBS) @@ -876,6 +892,7 @@ distclean-compile: @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/pkregann-pkregann.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/pksetmask.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/pksieve.Po@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/pkstat.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/pkstatascii.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/pkstatogr.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/pksvm.Po@am__quote@ @@ -1120,7 +1137,8 @@ maintainer-clean-generic: @echo "it deletes files that may require special tools to rebuild." clean: clean-am -clean-am: clean-binPROGRAMS clean-generic clean-libtool mostlyclean-am +clean-am: clean-binPROGRAMS clean-generic clean-libtool \ + clean-noinstPROGRAMS mostlyclean-am distclean: distclean-am -rm -rf ./$(DEPDIR) @@ -1191,15 +1209,16 @@ uninstall-am: uninstall-binPROGRAMS .MAKE: install-am install-strip .PHONY: CTAGS GTAGS TAGS all all-am check check-am clean \ - clean-binPROGRAMS clean-generic clean-libtool cscopelist-am \ - ctags ctags-am distclean distclean-compile distclean-generic \ - distclean-libtool distclean-tags distdir dvi dvi-am html \ - html-am info info-am install install-am install-binPROGRAMS \ - install-data install-data-am install-dvi install-dvi-am \ - install-exec install-exec-am install-html install-html-am \ - install-info install-info-am install-man install-pdf \ - install-pdf-am install-ps install-ps-am install-strip \ - installcheck installcheck-am installdirs maintainer-clean \ + clean-binPROGRAMS clean-generic clean-libtool \ + clean-noinstPROGRAMS cscopelist-am ctags ctags-am distclean \ + distclean-compile distclean-generic distclean-libtool \ + distclean-tags distdir dvi dvi-am html html-am info info-am \ + install install-am install-binPROGRAMS install-data \ + install-data-am install-dvi install-dvi-am install-exec \ + install-exec-am install-html install-html-am install-info \ + install-info-am install-man install-pdf install-pdf-am \ + install-ps install-ps-am install-strip installcheck \ + installcheck-am installdirs maintainer-clean \ maintainer-clean-generic mostlyclean mostlyclean-compile \ mostlyclean-generic mostlyclean-libtool pdf pdf-am ps ps-am \ tags tags-am uninstall uninstall-am uninstall-binPROGRAMS diff --git a/src/apps/pkann.cc b/src/apps/pkann.cc index 95cc3fb..932e3ce 100644 --- a/src/apps/pkann.cc +++ b/src/apps/pkann.cc @@ -817,51 +817,49 @@ int main(int argc, char *argv[]) } oldRowMask=rowMask; } - } - else - continue;//no coverage in this mask - short theMask=0; - for(short ivalue=0;ivalue<msknodata_opt.size();++ivalue){ - if(msknodata_opt[ivalue]>=0){//values set in msknodata_opt are invalid - if(lineMask[colMask]==msknodata_opt[ivalue]){ - theMask=lineMask[colMask]; - masked=true; - break; + short theMask=0; + for(short ivalue=0;ivalue<msknodata_opt.size();++ivalue){ + if(msknodata_opt[ivalue]>=0){//values set in msknodata_opt are invalid + if(lineMask[colMask]==msknodata_opt[ivalue]){ + theMask=lineMask[colMask]; + masked=true; + break; + } } - } - else{//only values set in msknodata_opt are valid - if(lineMask[colMask]!=-msknodata_opt[ivalue]){ - theMask=lineMask[colMask]; - masked=true; - } - else{ - masked=false; - break; + else{//only values set in msknodata_opt are valid + if(lineMask[colMask]!=-msknodata_opt[ivalue]){ + theMask=lineMask[colMask]; + masked=true; + } + else{ + masked=false; + break; + } } } + if(masked){ + if(classBag_opt.size()) + for(int ibag=0;ibag<nbag;++ibag) + classBag[ibag][icol]=theMask; + classOut[icol]=theMask; + continue; + } } - if(masked){ + bool valid=false; + for(int iband=0;iband<nband;++iband){ + if(hpixel[icol][iband]){ + valid=true; + break; + } + } + if(!valid){ if(classBag_opt.size()) for(int ibag=0;ibag<nbag;++ibag) - classBag[ibag][icol]=theMask; - classOut[icol]=theMask; - continue; + classBag[ibag][icol]=nodata_opt[0]; + classOut[icol]=nodata_opt[0]; + continue;//next column } } - bool valid=false; - for(int iband=0;iband<nband;++iband){ - if(hpixel[icol][iband]){ - valid=true; - break; - } - } - if(!valid){ - if(classBag_opt.size()) - for(int ibag=0;ibag<nbag;++ibag) - classBag[ibag][icol]=nodata_opt[0]; - classOut[icol]=nodata_opt[0]; - continue;//next column - } for(int iclass=0;iclass<nclass;++iclass) probOut[iclass][icol]=0; if(verbose_opt[0]>1) diff --git a/src/apps/pkcomposite.cc b/src/apps/pkcomposite.cc index 2c4aeb0..8fa7456 100644 --- a/src/apps/pkcomposite.cc +++ b/src/apps/pkcomposite.cc @@ -48,7 +48,7 @@ int main(int argc, char *argv[]) Optionpk<double> lry_opt("lry", "lry", "Lower right y value bounding box", 0.0); Optionpk<string> crule_opt("cr", "crule", "Composite rule (overwrite, maxndvi, maxband, minband, mean, mode (only for byte images), median, sum, maxallbands, minallbands", "overwrite"); Optionpk<int> ruleBand_opt("cb", "cband", "band index used for the composite rule (e.g., for ndvi, use --cband=0 --cband=1 with 0 and 1 indices for red and nir band respectively", 0); - Optionpk<double> srcnodata_opt("srcnodata", "srcnodata", "invalid value for input image", 0); + Optionpk<double> srcnodata_opt("srcnodata", "srcnodata", "invalid value for input image"); Optionpk<int> bndnodata_opt("bndnodata", "bndnodata", "Bands in input image to check if pixel is valid (used for srcnodata, min and max options)", 0); Optionpk<double> minValue_opt("min", "min", "flag values smaller or equal to this value as invalid."); Optionpk<double> maxValue_opt("max", "max", "flag values larger or equal to this value as invalid."); @@ -58,7 +58,7 @@ int main(int argc, char *argv[]) Optionpk<string> oformat_opt("of", "oformat", "Output image format (see also gdal_translate). Empty string: inherit from input image"); Optionpk<string> option_opt("co", "co", "Creation option for output file. Multiple options can be specified."); Optionpk<string> projection_opt("a_srs", "a_srs", "Override the spatial reference for the output file (leave blank to copy from input file, use epsg:3035 to use European projection and force to European grid"); - Optionpk<bool> file_opt("file", "file", "write number of observations for each pixels as additional layer in composite", false); + Optionpk<short> file_opt("file", "file", "write number of observations (1) or sequence nr of selected file (2) for each pixels as additional layer in composite", 0); Optionpk<short> weight_opt("w", "weight", "Weights (type: short) for the composite, use one weight for each input file in same order as input files are provided). Use value 1 for equal weights.", 1); Optionpk<short> class_opt("c", "class", "classes for multi-band output image: each band represents the number of observations for one specific class. Use value 0 for no multi-band output image.", 0); Optionpk<string> colorTable_opt("ct", "ct", "color table file with 5 columns: id R G B ALFA (0: transparent, 255: solid)"); @@ -131,8 +131,10 @@ int main(int argc, char *argv[]) cruleMap["maxallbands"]=crule::maxallbands; cruleMap["minallbands"]=crule::minallbands; - while(srcnodata_opt.size()<bndnodata_opt.size()) - srcnodata_opt.push_back(srcnodata_opt[0]); + if(srcnodata_opt.size()){ + while(srcnodata_opt.size()<bndnodata_opt.size()) + srcnodata_opt.push_back(srcnodata_opt[0]); + } while(bndnodata_opt.size()<srcnodata_opt.size()) bndnodata_opt.push_back(bndnodata_opt[0]); if(minValue_opt.size()){ @@ -614,6 +616,8 @@ int main(int argc, char *argv[]) break; } if(readValid){ + if(file_opt[0]==1) + ++fileBuffer[ib]; if(writeValid[ib]){ int iband=0; switch(cruleMap[crule_opt[0]]){ @@ -645,8 +649,8 @@ int main(int argc, char *argv[]) val_new=(readCol-0.5-lowerCol)*readBuffer[iband][upperCol-startCol]+(1-readCol+0.5+lowerCol)*readBuffer[iband][lowerCol-startCol]; writeBuffer[iband][ib]=val_new; } - fileBuffer[ib]=ifile; - // ++fileBuffer[ib]; + if(file_opt[0]>1) + fileBuffer[ib]=ifile; } break; default: @@ -660,8 +664,8 @@ int main(int argc, char *argv[]) val_new=readBuffer[iband][readCol-startCol]; writeBuffer[iband][ib]=val_new; } - fileBuffer[ib]=ifile; - // ++fileBuffer[ib]; + if(file_opt[0]>1) + fileBuffer[ib]=ifile; } break; } @@ -689,8 +693,8 @@ int main(int argc, char *argv[]) val_new*=weight_opt[ifile]; writeBuffer[iband][ib]=val_new; } - // fileBuffer[ib]=ifile; - ++fileBuffer[ib]; + if(file_opt[0]>1) + fileBuffer[ib]=ifile; } break; default: @@ -703,8 +707,8 @@ int main(int argc, char *argv[]) val_new*=weight_opt[ifile]; writeBuffer[iband][ib]=val_new; } - // fileBuffer[ib]=ifile; - ++fileBuffer[ib]; + if(file_opt[0]>1) + fileBuffer[ib]=ifile; } break; } @@ -768,7 +772,8 @@ int main(int argc, char *argv[]) } break; } - ++fileBuffer[ib]; + if(file_opt[0]>1) + fileBuffer[ib]=ifile; break; case(crule::overwrite): default: @@ -797,8 +802,8 @@ int main(int argc, char *argv[]) } break; } - // fileBuffer[ib]=ifile; - ++fileBuffer[ib]; + if(file_opt[0]>1) + fileBuffer[ib]=ifile; break; } } @@ -836,8 +841,9 @@ int main(int argc, char *argv[]) } break; } - ++fileBuffer[ib]; - break; + if(file_opt[0]>1) + fileBuffer[ib]=ifile; + break; case(crule::mode): switch(theResample){ case(BILINEAR): @@ -891,8 +897,8 @@ int main(int argc, char *argv[]) } break; } - // fileBuffer[ib]=ifile; - ++fileBuffer[ib]; + if(file_opt[0]>1) + fileBuffer[ib]=ifile; break; } } @@ -920,7 +926,8 @@ int main(int argc, char *argv[]) vector<short>::iterator maxit=maxBuffer[icol].begin(); maxit=stat.mymax(maxBuffer[icol],maxBuffer[icol].begin(),maxBuffer[icol].end()); writeBuffer[0][icol]=distance(maxBuffer[icol].begin(),maxit); - fileBuffer[icol]=*(maxit); + if(file_opt[0]>1) + fileBuffer[icol]=*(maxit); } try{ imgWriter.writeData(writeBuffer[0],GDT_Float64,irow,0); @@ -940,27 +947,32 @@ int main(int argc, char *argv[]) for(int icol=0;icol<imgWriter.nrOfCol();++icol){ switch(cruleMap[crule_opt[0]]){ case(crule::mean): - assert(storeBuffer[bands[iband]][icol].size()==fileBuffer[icol]); + if(file_opt[0]<2) + // assert(storeBuffer[bands[iband]][icol].size()==fileBuffer[icol]); if(storeBuffer[bands[iband]][icol].size()) writeBuffer[iband][icol]=stat.mean(storeBuffer[bands[iband]][icol]); break; case(crule::median): - assert(storeBuffer[bands[iband]][icol].size()==fileBuffer[icol]); + if(file_opt[0]<2) + // assert(storeBuffer[bands[iband]][icol].size()==fileBuffer[icol]); if(storeBuffer[bands[iband]][icol].size()) writeBuffer[iband][icol]=stat.median(storeBuffer[bands[iband]][icol]); break; case(crule::sum)://sum - assert(storeBuffer[bands[iband]][icol].size()==fileBuffer[icol]); + if(file_opt[0]<2) + // assert(storeBuffer[bands[iband]][icol].size()==fileBuffer[icol]); if(storeBuffer[bands[iband]][icol].size()) writeBuffer[iband][icol]=stat.sum(storeBuffer[bands[iband]][icol]); break; case(crule::minallbands): - assert(storeBuffer[bands[iband]][icol].size()==fileBuffer[icol]); + if(file_opt[0]<2) + // assert(storeBuffer[bands[iband]][icol].size()==fileBuffer[icol]); if(storeBuffer[bands[iband]][icol].size()) writeBuffer[iband][icol]=stat.mymin(storeBuffer[bands[iband]][icol]); break; case(crule::maxallbands): - assert(storeBuffer[bands[iband]][icol].size()==fileBuffer[icol]); + if(file_opt[0]<2) + // assert(storeBuffer[bands[iband]][icol].size()==fileBuffer[icol]); if(storeBuffer[bands[iband]][icol].size()) writeBuffer[iband][icol]=stat.mymax(storeBuffer[bands[iband]][icol]); break; diff --git a/src/apps/pkcrop.cc b/src/apps/pkcrop.cc index 6d6c1e9..526a498 100644 --- a/src/apps/pkcrop.cc +++ b/src/apps/pkcrop.cc @@ -457,7 +457,12 @@ int main(int argc, char *argv[]) double theMin=0; double theMax=0; if(autoscale_opt.size()){ - imgReader.getMinMax(static_cast<int>(startCol),static_cast<int>(endCol),static_cast<int>(startRow),static_cast<int>(endRow),readBand,theMin,theMax); + try{ + imgReader.getMinMax(static_cast<int>(startCol),static_cast<int>(endCol),static_cast<int>(startRow),static_cast<int>(endRow),readBand,theMin,theMax); + } + catch(string errorString){ + cout << errorString << endl; + } if(verbose_opt[0]) cout << "minmax: " << theMin << ", " << theMax << endl; double theScale=(autoscale_opt[1]-autoscale_opt[0])/(theMax-theMin); diff --git a/src/apps/pkdiff.cc b/src/apps/pkdiff.cc index 375c591..39a2643 100644 --- a/src/apps/pkdiff.cc +++ b/src/apps/pkdiff.cc @@ -34,7 +34,7 @@ int main(int argc, char *argv[]) Optionpk<string> layer_opt("ln", "ln", "Layer name(s) in sample. Leave empty to select all (for vector reference datasets only)"); 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."); 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); - Optionpk<short> nodata_opt("nodata", "nodata", "No data value(s) in input or reference dataset are ignored"); + Optionpk<double> nodata_opt("nodata", "nodata", "No data value(s) in input or reference dataset are ignored"); 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); Optionpk<bool> rmse_opt("rmse", "rmse", "Report root mean squared error", false); Optionpk<bool> regression_opt("reg", "reg", "Report linear regression (Input = c0+c1*Reference)", false); @@ -168,7 +168,7 @@ int main(int argc, char *argv[]) for(int iflag=0;iflag<nodata_opt.size();++iflag){ vector<short>::iterator fit; - fit=find(inputRange.begin(),inputRange.end(),nodata_opt[iflag]); + fit=find(inputRange.begin(),inputRange.end(),static_cast<short>(nodata_opt[iflag])); if(fit!=inputRange.end()) inputRange.erase(fit); } @@ -344,13 +344,13 @@ int main(int argc, char *argv[]) } x=poPoint->getX(); y=poPoint->getY(); - short inputValue; - vector<short> inputValues; + double inputValue; + vector<double> inputValues; bool isHomogeneous=true; short maskValue; short outputValue; //read referenceValue from feature - short referenceValue; + unsigned short referenceValue; string referenceClassName; if(classValueMap.size()){ referenceClassName=readFeature->GetFieldAsString(readFeature->GetFieldIndex(labelref_opt[0].c_str())); @@ -476,12 +476,12 @@ int main(int argc, char *argv[]) ++ntotalValidation; if(classValueMap.size()){ assert(inputValue<nameVector.size()); - string className=nameVector[inputValue]; + string className=nameVector[static_cast<unsigned short>(inputValue)]; cm.incrementResult(type2string<short>(classValueMap[referenceClassName]),type2string<short>(classValueMap[className]),1); } else{ - int rc=distance(referenceRange.begin(),find(referenceRange.begin(),referenceRange.end(),referenceValue)); - int ic=distance(inputRange.begin(),find(inputRange.begin(),inputRange.end(),inputValue)); + int rc=distance(referenceRange.begin(),find(referenceRange.begin(),referenceRange.end(),static_cast<unsigned short>(referenceValue))); + int ic=distance(inputRange.begin(),find(inputRange.begin(),inputRange.end(),static_cast<unsigned short>(inputValue))); assert(rc<nclass); assert(ic<nclass); ++nvalidation[rc]; @@ -524,18 +524,18 @@ int main(int argc, char *argv[]) else fs << labelclass_opt[0]; if(output_opt.size()) - writeFeature->SetField(fs.str().c_str(),static_cast<int>(inputValue)); + writeFeature->SetField(fs.str().c_str(),inputValue); if(!windowJ&&!windowI){//centre pixel if(confusion_opt[0]){ ++ntotalValidation; if(classValueMap.size()){ assert(inputValue<nameVector.size()); - string className=nameVector[inputValue]; + string className=nameVector[static_cast<unsigned short>(inputValue)]; cm.incrementResult(type2string<short>(classValueMap[referenceClassName]),type2string<short>(classValueMap[className]),1); } else{ - int rc=distance(referenceRange.begin(),find(referenceRange.begin(),referenceRange.end(),referenceValue)); - int ic=distance(inputRange.begin(),find(inputRange.begin(),inputRange.end(),inputValue)); + int rc=distance(referenceRange.begin(),find(referenceRange.begin(),referenceRange.end(),static_cast<unsigned short>(referenceValue))); + int ic=distance(inputRange.begin(),find(inputRange.begin(),inputRange.end(),static_cast<unsigned short>(inputValue))); if(rc>=nclass) continue; if(ic>=nclass) @@ -656,7 +656,7 @@ int main(int argc, char *argv[]) referenceReaderGdal.getRange(referenceRange,band_opt[1]); for(int iflag=0;iflag<nodata_opt.size();++iflag){ vector<short>::iterator fit; - fit=find(referenceRange.begin(),referenceRange.end(),nodata_opt[iflag]); + fit=find(referenceRange.begin(),referenceRange.end(),static_cast<unsigned short>(nodata_opt[iflag])); if(fit!=referenceRange.end()) referenceRange.erase(fit); } @@ -856,7 +856,8 @@ int main(int argc, char *argv[]) }//raster dataset if(confusion_opt[0]){ - assert(cm.nReference()); + + // assert(cm.nReference()); cout << cm << endl; cout << "class #samples userAcc prodAcc" << endl; double se95_ua=0; diff --git a/src/apps/pkdumpimg.cc b/src/apps/pkdumpimg.cc index c6334e6..149f50f 100644 --- a/src/apps/pkdumpimg.cc +++ b/src/apps/pkdumpimg.cc @@ -166,8 +166,8 @@ int main(int argc, char *argv[]) extentReader.close(); } } - double uli,ulj,lri,lrj;//image coordinates - double magicX=1,magicY=1; + double uli,ulj,lri,lrj;//image coordinates + double magicX=1,magicY=1; if(ulx_opt[0]>=lrx_opt[0]){//default bounding box: no cropping uli=0; lri=imgReader.nrOfCol()-1; @@ -178,17 +178,15 @@ int main(int argc, char *argv[]) imgReader.getBoundingBox(cropulx,cropuly,croplrx,croplry); imgReader.geo2image(cropulx+(magicX-1.0)*imgReader.getDeltaX(),cropuly-(magicY-1.0)*imgReader.getDeltaY(),uli,ulj); imgReader.geo2image(croplrx+(magicX-2.0)*imgReader.getDeltaX(),croplry-(magicY-2.0)*imgReader.getDeltaY(),lri,lrj); - ncropcol=ceil((croplrx-cropulx)/dx); - ncroprow=ceil((cropuly-croplry)/dy); - ncropcol=ceil((croplrx-cropulx)/dx); - ncroprow=ceil((cropuly-croplry)/dy); + ncropcol=abs(static_cast<int>(ceil((croplrx-cropulx)/dx))); + ncroprow=abs(static_cast<int>(ceil((cropuly-croplry)/dy))); } else{ imgReader.geo2image(cropulx+(magicX-1.0)*imgReader.getDeltaX(),cropuly-(magicY-1.0)*imgReader.getDeltaY(),uli,ulj); imgReader.geo2image(croplrx+(magicX-2.0)*imgReader.getDeltaX(),croplry-(magicY-2.0)*imgReader.getDeltaY(),lri,lrj); - ncropcol=ceil((croplrx-cropulx)/dx); - ncroprow=ceil((cropuly-croplry)/dy); + ncropcol=abs(static_cast<int>(ceil((croplrx-cropulx)/dx))); + ncroprow=abs(static_cast<int>(ceil((cropuly-croplry)/dy))); uli=floor(uli); ulj=floor(ulj); lri=floor(lri); diff --git a/src/apps/pkenhance.cc b/src/apps/pkenhance.cc index 7414156..200ae1d 100644 --- a/src/apps/pkenhance.cc +++ b/src/apps/pkenhance.cc @@ -138,8 +138,8 @@ int main(int argc,char **argv) { //calculate histograms unsigned int nbinRef=nbin_opt[0]; unsigned int nbinInput=nbin_opt[0]; - std::vector<unsigned long int> histRef(nbinRef); - std::vector<unsigned long int> histInput(nbinInput); + std::vector<double> histRef(nbinRef); + std::vector<double> histInput(nbinInput); double minValueRef=0; double maxValueRef=0; double minValueInput=0; diff --git a/src/apps/pkextract.cc b/src/apps/pkextract.cc index 6aea515..de7eabe 100644 --- a/src/apps/pkextract.cc +++ b/src/apps/pkextract.cc @@ -1232,10 +1232,6 @@ int main(int argc, char *argv[]) double theValue=0; for(int index=0;index<windowValues.size();++index){ try{ - if(windowValues[index].size()!=buffer_opt[0]*buffer_opt[0]){ - std::string errorString="windowValues[index].size()!=buffer*buffer"; - throw(errorString); - } if(ruleMap[rule_opt[0]]==rule::mean) theValue=stat.mean(windowValues[index]); else if(ruleMap[rule_opt[0]]==rule::stdev) diff --git a/src/apps/pkinfo.cc b/src/apps/pkinfo.cc index e244675..9e433cb 100644 --- a/src/apps/pkinfo.cc +++ b/src/apps/pkinfo.cc @@ -43,9 +43,6 @@ int main(int argc, char *argv[]) Optionpk<bool> min_opt("min", "minimum", "Shows min value of the image ", false,0); Optionpk<bool> max_opt("max", "maximum", "Shows max value of the image ", false,0); Optionpk<bool> stat_opt("stats", "statistics", "Shows statistics (min,max, mean and stdDev of the image)", false,0); - Optionpk<double> src_min_opt("src_min", "src_min", "Sets minimum for histogram (does not calculate min value: use -mm instead)"); - Optionpk<double> src_max_opt("src_max", "src_max", "Sets maximum for histogram (does not calculate min value: use -mm instead)"); - Optionpk<bool> relative_opt("rel", "relative", "Calculates relative histogram in percentage", false,0); Optionpk<bool> projection_opt("a_srs", "a_srs", "Shows projection of the image ", false,0); Optionpk<bool> geo_opt("geo", "geo", "Gets geotransform ", false,0); Optionpk<bool> interleave_opt("il", "interleave", "Shows interleave ", false,0); @@ -61,8 +58,6 @@ int main(int argc, char *argv[]) Optionpk<double> uly_opt("uly", "uly", "Upper left y value bounding box"); Optionpk<double> lrx_opt("lrx", "lrx", "Lower right x value bounding box"); Optionpk<double> lry_opt("lry", "lry", "Lower right y value bounding box"); - Optionpk<bool> hist_opt("hist", "histogram", "Calculates histogram. Use --rel for a relative histogram output. ", false,0); - Optionpk<unsigned int> nbin_opt("nbin", "nbin", "Number of bins used in histogram. Use 0 for all input values as integers"); Optionpk<bool> type_opt("ot", "otype", "Returns data type", false,0); Optionpk<bool> description_opt("d", "description", "Returns image description", false,0); Optionpk<bool> metadata_opt("meta", "meta", "Shows meta data ", false,0); @@ -85,9 +80,6 @@ int main(int argc, char *argv[]) min_opt.retrieveOption(argc,argv); max_opt.retrieveOption(argc,argv); stat_opt.retrieveOption(argc,argv); - src_min_opt.retrieveOption(argc,argv); - src_max_opt.retrieveOption(argc,argv); - relative_opt.retrieveOption(argc,argv); projection_opt.retrieveOption(argc,argv); geo_opt.retrieveOption(argc,argv); interleave_opt.retrieveOption(argc,argv); @@ -103,8 +95,6 @@ int main(int argc, char *argv[]) uly_opt.retrieveOption(argc,argv); lrx_opt.retrieveOption(argc,argv); lry_opt.retrieveOption(argc,argv); - hist_opt.retrieveOption(argc,argv); - nbin_opt.retrieveOption(argc,argv); type_opt.retrieveOption(argc,argv); description_opt.retrieveOption(argc,argv); metadata_opt.retrieveOption(argc,argv); @@ -296,47 +286,6 @@ int main(int argc, char *argv[]) } } } - if(relative_opt[0]) - hist_opt[0]=true; - if(hist_opt[0]){ - assert(band_opt[0]<imgReader.nrOfBand()); - unsigned int nbin=(nbin_opt.size())? nbin_opt[0]:0; - std::vector<unsigned long int> output; - minValue=0; - maxValue=0; - //todo: optimize such that getMinMax is only called once... - imgReader.getMinMax(minValue,maxValue,band_opt[0]); - - if(src_min_opt.size()) - minValue=src_min_opt[0]; - if(src_max_opt.size()) - maxValue=src_max_opt[0]; - unsigned long int nsample=imgReader.getHistogram(output,minValue,maxValue,nbin,band_opt[0]); - std::cout.precision(10); - for(int bin=0;bin<nbin;++bin){ - double binValue=0; - if(nbin==maxValue-minValue+1) - binValue=minValue+bin; - else - binValue=minValue+static_cast<double>(maxValue-minValue)*(bin+0.5)/nbin; - std::cout << binValue << " "; - if(relative_opt[0]) - std::cout << 100.0*static_cast<double>(output[bin])/static_cast<double>(nsample) << std::endl; - else - std::cout << static_cast<double>(output[bin]) << std::endl; - } - } - // else{ - // int minCol,minRow; - // if(src_min_opt.size()){ - // assert(band_opt[0]<imgReader.nrOfBand()); - // std::cout << "-min " << imgReader.getMin(minCol, minRow,band_opt[0]); - // } - // if(src_max_opt.size()){ - // assert(band_opt[0]<imgReader.nrOfBand()); - // std::cout << "-max " << imgReader.getMax(minCol, minRow,band_opt[0]); - // } - // } if(projection_opt[0]){ if(imgReader.isGeoRef()) std::cout << " -a_srs " << imgReader.getProjection() << " "; @@ -429,6 +378,6 @@ int main(int argc, char *argv[]) else std::cout << "no intersect" << std::endl; } - if(!read_opt[0]&&!hist_opt[0]) + if(!read_opt[0]) std::cout << std::endl; } diff --git a/src/apps/pkkalman.cc b/src/apps/pkkalman.cc index 96885db..3d5cc7c 100644 --- a/src/apps/pkkalman.cc +++ b/src/apps/pkkalman.cc @@ -43,23 +43,27 @@ int main(int argc,char **argv) { Optionpk<string> outputfb_opt("ofb", "outputfb", "Output raster dataset for smooth model"); Optionpk<double> modnodata_opt("modnodata", "modnodata", "invalid value for model input", 0); Optionpk<double> obsnodata_opt("obsnodata", "obsnodata", "invalid value for observation input", 0); - Optionpk<double> modoffset_opt("modoffset", "modoffset", "offset used to read model input dataset (value=offset+scale*readValue", 0); - Optionpk<double> obsoffset_opt("obsoffset", "obsoffset", "offset used to read observation input dataset (value=offset+scale*readValue", 0); - Optionpk<double> modscale_opt("modscale", "modscale", "scale used to read model input dataset (value=offset+scale*readValue", 1); - Optionpk<double> obsscale_opt("obsscale", "obsscale", "scale used to read observation input dataset (value=offset+scale*readValue", 1); + Optionpk<double> modoffset_opt("modoffset", "modoffset", "offset used to read model input dataset (value=offset+scale*readValue"); + Optionpk<double> obsoffset_opt("obsoffset", "obsoffset", "offset used to read observation input dataset (value=offset+scale*readValue"); + Optionpk<double> modscale_opt("modscale", "modscale", "scale used to read model input dataset (value=offset+scale*readValue"); + Optionpk<double> obsscale_opt("obsscale", "obsscale", "scale used to read observation input dataset (value=offset+scale*readValue"); Optionpk<double> eps_opt("eps", "eps", "epsilon for non zero division", 0.00001); Optionpk<double> uncertModel_opt("um", "uncertmodel", "Multiply this value with std dev of first model image to obtain uncertainty of model",2); Optionpk<double> uncertObs_opt("uo", "uncertobs", "Uncertainty of valid observations",0); + Optionpk<double> weight_opt("w", "weight", "Set observation uncertainty as weighted difference between observation and model (use -w 0 to use a constant observation uncertainty, use -w value >> 1 to penalize low observation values with respect to model, use -w value << 0 to penalize a high observation values with respect to model"); Optionpk<double> deltaObs_opt("dobs", "deltaobs", "Lower and upper thresholds for relative pixel differences (in percentage): (observation-model)/model. For instance to force the observation within a +/- 10 % interval, use: -dobs -10 -dobs 10 (equivalent to -dobs 10). Leave empty to always update on observation"); Optionpk<double> uncertNodata_opt("unodata", "uncertnodata", "Uncertainty in case of no-data values in observation", 10000); // Optionpk<double> regTime_opt("rt", "regtime", "Relative Weight for regression in time series", 1.0); Optionpk<bool> regSensor_opt("rs", "regsensor", "Set optional regression for sensor difference (model - observation).", false); - Optionpk<int> down_opt("down", "down", "Downsampling factor for reading model data to calculate regression", 9); + Optionpk<int> down_opt("down", "down", "Downsampling factor for reading model data to calculate regression"); Optionpk<float> threshold_opt("th", "threshold", "threshold for selecting samples (randomly). Provide probability in percentage (>0) or absolute (<0).", 0); Optionpk<int> minreg_opt("minreg", "minreg", "Minimum number of pixels to take into account for regression", 5, 2); // Optionpk<bool> regObs_opt("regobs", "regobs", "Perform regression between modeled and observed value",false); // Optionpk<double> checkDiff_opt("diff", "diff", "Flag observation as invalid if difference with model is above uncertainty",false); Optionpk<unsigned short> window_opt("win", "window", "window size for calculating regression (use 0 for global)", 0); + // 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."); + // Optionpk<double> msknodata_opt("msknodata", "msknodata", "Mask value(s) not to consider for filtering. First value will be set in output image.", 0); + Optionpk<string> oformat_opt("of", "oformat", "Output image format (see also gdal_translate). Empty string: inherit from input image","GTiff",2); Optionpk<string> option_opt("co", "co", "Creation option for output file. Multiple options can be specified."); Optionpk<short> verbose_opt("v", "verbose", "verbose mode when positive", 0); @@ -84,6 +88,7 @@ int main(int argc,char **argv) { eps_opt.retrieveOption(argc,argv); uncertModel_opt.retrieveOption(argc,argv); uncertObs_opt.retrieveOption(argc,argv); + weight_opt.retrieveOption(argc,argv); deltaObs_opt.retrieveOption(argc,argv); uncertNodata_opt.retrieveOption(argc,argv); // regTime_opt.retrieveOption(argc,argv); @@ -94,6 +99,8 @@ int main(int argc,char **argv) { // regObs_opt.retrieveOption(argc,argv); // checkDiff_opt.retrieveOption(argc,argv); window_opt.retrieveOption(argc,argv); + // mask_opt.retrieveOption(argc,argv); + // msknodata_opt.retrieveOption(argc,argv); oformat_opt.retrieveOption(argc,argv); option_opt.retrieveOption(argc,argv); verbose_opt.retrieveOption(argc,argv); @@ -107,6 +114,13 @@ int main(int argc,char **argv) { exit(0);//help was invoked, stop processing } + if(deltaObs_opt.size()==1){ + if(deltaObs_opt[0]<=0) + deltaObs_opt.push_back(-deltaObs_opt[0]); + else + deltaObs_opt.insert(deltaObs_opt.begin(),-deltaObs_opt[0]); + } + try{ ostringstream errorStream; if(model_opt.size()<2){ @@ -117,12 +131,6 @@ int main(int argc,char **argv) { errorStream << "Error: no observation dataset selected, use option -obs" << endl; throw(errorStream.str()); } - if(deltaObs_opt.size()==1){ - if(deltaObs_opt[0]<=0) - deltaObs_opt.push_back(-deltaObs_opt[0]); - else - deltaObs_opt.insert(deltaObs_opt.begin(),-deltaObs_opt[0]); - } if(direction_opt[0]=="smooth"){ if(outputfw_opt.empty()){ errorStream << "Error: output forward datasets is not provided, use option -ofw" << endl; @@ -199,8 +207,39 @@ int main(int argc,char **argv) { option_opt.push_back(theInterleave); } + if(down_opt.empty()){ + imgReaderModel1.open(model_opt[0]); + double resModel=imgReaderModel1.getDeltaX(); + double resObs=imgReaderObs.getDeltaX(); + int down=static_cast<int>(ceil(resModel/resObs)); + if(!(down%2)) + down+=1; + down_opt.push_back(down); + imgReaderModel1.close(); + } imgReaderObs.close(); + // ImgReaderGdal maskReader; + // double colMask=0; + // double rowMask=0; + + // if(mask_opt.size()){ + // try{ + // if(verbose_opt[0]>=1) + // std::cout << "opening mask image file " << mask_opt[0] << std::endl; + // maskReader.open(mask_opt[0]); + // maskReader.setNoData(msknodata_opt); + // } + // catch(string error){ + // cerr << error << std::endl; + // exit(2); + // } + // catch(...){ + // cerr << "error catched" << std::endl; + // exit(1); + // } + // } + int obsindex=0; const char* pszMessage; @@ -237,6 +276,9 @@ int main(int argc,char **argv) { // cout << "tobservation_opt[tindex] " << tobservation_opt[tindex] << " " << relobsindex.back() << endl; } + double geox=0; + double geoy=0; + if(find(direction_opt.begin(),direction_opt.end(),"forward")!=direction_opt.end()){ ///////////////////////////// forward model ///////////////////////// cout << "Running forward model" << endl; @@ -268,8 +310,10 @@ int main(int argc,char **argv) { try{ imgReaderModel1.open(model_opt[0]); imgReaderModel1.setNoData(modnodata_opt); - imgReaderModel1.setOffset(modoffset_opt[0]); - imgReaderModel1.setScale(modscale_opt[0]); + if(modoffset_opt.size()) + imgReaderModel1.setOffset(modoffset_opt[0]); + if(modscale_opt.size()) + imgReaderModel1.setScale(modscale_opt[0]); } catch(string errorString){ cerr << errorString << endl; @@ -284,8 +328,6 @@ int main(int argc,char **argv) { double minValue, maxValue, meanValue, stdDev; void* pProgressData; rasterBand->ComputeStatistics(0,&minValue,&maxValue,&meanValue,&stdDev,pfnProgress,pProgressData); - double x=0; - double y=0; double modRow=0; double modCol=0; if(relobsindex[0]>0){//initialize output_opt[0] as model[0] @@ -296,17 +338,46 @@ int main(int argc,char **argv) { vector<double> estReadBuffer; vector<double> estWriteBuffer(ncol); vector<double> uncertWriteBuffer(ncol); - imgWriterEst.image2geo(0,irow,x,y); - imgReaderModel1.geo2image(x,y,modCol,modRow); + // vector<double> lineMask; + imgWriterEst.image2geo(0,irow,geox,geoy); + imgReaderModel1.geo2image(geox,geoy,modCol,modRow); assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow()); try{ imgReaderModel1.readData(estReadBuffer,GDT_Float64,modRow); //simple nearest neighbor //stat.nearUp(estReadBuffer,estWriteBuffer); + // double oldRowMask=-1;//keep track of row mask to optimize number of line readings for(int icol=0;icol<ncol;++icol){ - imgWriterEst.image2geo(icol,irow,x,y); - imgReaderModel1.geo2image(x,y,modCol,modRow); + imgWriterEst.image2geo(icol,irow,geox,geoy); + // bool masked=false; + // if(mask_opt.size()){ + // //read mask + // maskReader.geo2image(geox,geoy,colMask,rowMask); + // colMask=static_cast<int>(colMask); + // rowMask=static_cast<int>(rowMask); + // if(rowMask>=0&&rowMask<maskReader.nrOfRow()&&colMask>=0&&colMask<maskReader.nrOfCol()){ + // if(static_cast<int>(rowMask)!=static_cast<int>(oldRowMask)){ + // try{ + // maskReader.readData(lineMask,GDT_Int16,static_cast<int>(rowMask)); + // } + // catch(string errorstring){ + // cerr << errorstring << endl; + // exit(1); + // } + // catch(...){ + // cerr << "error catched" << std::endl; + // exit(3); + // } + // oldRowMask=rowMask; + // } + // masked=(maskReader.isNoData(lineMask[colMask])); + // } + // estWriteBuffer[icol]=msknodata_opt[0]; + // uncertWriteBuffer[icol]=msknodata_opt[0]; + // continue;//next column + // } + imgReaderModel1.geo2image(geox,geoy,modCol,modRow); double modValue=estReadBuffer[modCol]; if(imgReaderModel1.isNoData(modValue)){ estWriteBuffer[icol]=obsnodata_opt[0]; @@ -334,35 +405,68 @@ int main(int argc,char **argv) { imgReaderObs.open(observation_opt[0]); imgReaderObs.getGeoTransform(geotransform); imgReaderObs.setNoData(obsnodata_opt); - imgReaderObs.setOffset(obsoffset_opt[0]); - imgReaderObs.setScale(obsscale_opt[0]); + if(obsoffset_opt.size()) + imgReaderObs.setOffset(obsoffset_opt[0]); + if(obsscale_opt.size()) + imgReaderObs.setScale(obsscale_opt[0]); if(regSensor_opt[0]) - errObs=imgreg.getRMSE(imgReaderModel1,imgReaderObs,c0obs,c1obs,verbose_opt[0]); + errObs=imgreg.getRMSE(imgReaderModel1,imgReaderObs,c0obs,c1obs,0,0,verbose_opt[0]); else{ c0obs=0; c1obs=1; errObs=0; } + if(verbose_opt[0]) + cout << "c0obs, c1obs: " << c0obs << ", " << c1obs << endl; for(int irow=0;irow<nrow;++irow){ vector<double> estReadBuffer; - imgWriterEst.image2geo(0,irow,x,y); - imgReaderModel1.geo2image(x,y,modCol,modRow); + imgWriterEst.image2geo(0,irow,geox,geoy); + imgReaderModel1.geo2image(geox,geoy,modCol,modRow); assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow()); imgReaderModel1.readData(estReadBuffer,GDT_Float64,modRow); vector<double> obsLineBuffer; vector<double> estWriteBuffer(ncol); vector<double> uncertWriteBuffer(ncol); vector<double> uncertObsLineBuffer; + // vector<double> lineMask; // imgReaderObs.readData(estWriteBuffer,GDT_Float64,irow,0); imgReaderObs.readData(obsLineBuffer,GDT_Float64,irow,0); if(imgReaderObs.nrOfBand()>1) imgReaderObs.readData(uncertObsLineBuffer,GDT_Float64,irow,1); + // double oldRowMask=-1;//keep track of row mask to optimize number of line readings for(int icol=0;icol<ncol;++icol){ - imgWriterEst.image2geo(icol,irow,x,y); - imgReaderModel1.geo2image(x,y,modCol,modRow); + imgWriterEst.image2geo(icol,irow,geox,geoy); + // bool masked=false; + // if(mask_opt.size()){ + // //read mask + // maskReader.geo2image(geox,geoy,colMask,rowMask); + // colMask=static_cast<int>(colMask); + // rowMask=static_cast<int>(rowMask); + // if(rowMask>=0&&rowMask<maskReader.nrOfRow()&&colMask>=0&&colMask<maskReader.nrOfCol()){ + // if(static_cast<int>(rowMask)!=static_cast<int>(oldRowMask)){ + // try{ + // maskReader.readData(lineMask,GDT_Int16,static_cast<int>(rowMask)); + // } + // catch(string errorstring){ + // cerr << errorstring << endl; + // exit(1); + // } + // catch(...){ + // cerr << "error catched" << std::endl; + // exit(3); + // } + // oldRowMask=rowMask; + // } + // masked=(maskReader.isNoData(lineMask[colMask])); + // } + // estWriteBuffer[icol]=msknodata_opt[0]; + // uncertWriteBuffer[icol]=msknodata_opt[0]; + // continue;//next column + // } + imgReaderModel1.geo2image(geox,geoy,modCol,modRow); assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow()); double modValue=estReadBuffer[modCol]; if(imgReaderModel1.isNoData(modValue)){//model is nodata: retain observation @@ -402,8 +506,9 @@ int main(int argc,char **argv) { if(!imgReaderObs.isNoData(obsLineBuffer[icol])){ double kalmanGain=1; double uncertObs=uncertObs_opt[0]; - bool doUpdate=true; - if(deltaObs_opt.size()){ + if(uncertObsLineBuffer.size()>icol) + uncertObs=uncertObsLineBuffer[icol]; + else if(weight_opt.size()||deltaObs_opt.size()){ vector<double> obsWindowBuffer;//buffer for observation to calculate average corresponding to model pixel int minCol=(icol>down_opt[0]/2) ? icol-down_opt[0]/2 : 0; int maxCol=(icol+down_opt[0]/2<imgReaderObs.nrOfCol()) ? icol+down_opt[0]/2 : imgReaderObs.nrOfCol()-1; @@ -415,29 +520,26 @@ int main(int argc,char **argv) { statobs.setNoDataValues(obsnodata_opt); double obsMeanValue=(statobs.mean(obsWindowBuffer)-c0obs)/c1obs; double difference=obsMeanValue-modValue; - if(modValue){ - difference/=modValue;//relative difference - difference*=100;//in percentage - doUpdate=(difference>=deltaObs_opt[0]);//lower bound - doUpdate&=(difference<=deltaObs_opt[1]);//upper bound - //todo: use deltaObs to set observation uncertainty instead of setting doUpdate - if(-difference>uncertObs_opt[0]) - uncertObs=-difference; - } - if(verbose_opt[0]>1){ - if(!doUpdate) - cout << "obsMeanValue:" << obsMeanValue << ", modValue: " << modValue << endl; - } + if(modValue&&deltaObs_opt.size()){ + double relativeDifference=100.0*difference/modValue; + if(relativeDifference<deltaObs_opt[0])//lower bound + kalmanGain=0; + else if(relativeDifference>deltaObs_opt[1])//upper bound + kalmanGain=0; + } + uncertObs=-weight_opt[0]*difference; + if(uncertObs<=0) + uncertObs=0; + if(verbose_opt[0]>1) + cout << "obsMeanValue:" << obsMeanValue << ", modValue: " << modValue << endl; } - if(doUpdate){ - if(uncertObsLineBuffer.size()>icol) - uncertObs=uncertObsLineBuffer[icol]; + if(kalmanGain>0){ if((uncertWriteBuffer[icol]+uncertObs)>eps_opt[0]) kalmanGain=uncertWriteBuffer[icol]/(uncertWriteBuffer[icol]+uncertObs); - assert(kalmanGain<=1); - estWriteBuffer[icol]+=kalmanGain*(obsLineBuffer[icol]-estWriteBuffer[icol]); - uncertWriteBuffer[icol]*=(1-kalmanGain); } + assert(kalmanGain<=1); + estWriteBuffer[icol]+=kalmanGain*(obsLineBuffer[icol]-estWriteBuffer[icol]); + uncertWriteBuffer[icol]*=(1-kalmanGain); } } imgWriterEst.writeData(estWriteBuffer,GDT_Float64,irow,0); @@ -476,12 +578,16 @@ int main(int argc,char **argv) { //calculate regression between two subsequence model inputs imgReaderModel1.open(model_opt[modindex-1]); imgReaderModel1.setNoData(modnodata_opt); - imgReaderModel1.setOffset(modoffset_opt[0]); - imgReaderModel1.setScale(modscale_opt[0]); + if(modoffset_opt.size()) + imgReaderModel1.setOffset(modoffset_opt[0]); + if(modscale_opt.size()) + imgReaderModel1.setScale(modscale_opt[0]); imgReaderModel2.open(model_opt[modindex]); imgReaderModel2.setNoData(modnodata_opt); - imgReaderModel2.setOffset(modoffset_opt[0]); - imgReaderModel2.setScale(modscale_opt[0]); + if(modoffset_opt.size()) + imgReaderModel2.setOffset(modoffset_opt[0]); + if(modscale_opt.size()) + imgReaderModel2.setScale(modscale_opt[0]); //calculate regression //we could re-use the points from second image from last run, but //to keep it general, we must redo it (overlap might have changed) @@ -491,8 +597,10 @@ int main(int argc,char **argv) { if(verbose_opt[0]) cout << "Calculating regression for " << imgReaderModel1.getFileName() << " " << imgReaderModel2.getFileName() << endl; - double errMod=imgreg.getRMSE(imgReaderModel1,imgReaderModel2,c0modGlobal,c1modGlobal); + double errMod=imgreg.getRMSE(imgReaderModel1,imgReaderModel2,c0modGlobal,c1modGlobal,0,0); // double errMod=imgreg.getRMSE(imgReaderModel1,imgReaderModel2,c0modGlobal,c1modGlobal,verbose_opt[0]); + if(verbose_opt[0]) + cout << "c0modGlobal, c1modGlobal: " << c0modGlobal << ", " << c1modGlobal << endl; bool update=false; if(obsindex<relobsindex.size()){ @@ -505,13 +613,15 @@ int main(int argc,char **argv) { imgReaderObs.open(observation_opt[obsindex]); imgReaderObs.getGeoTransform(geotransform); imgReaderObs.setNoData(obsnodata_opt); - imgReaderObs.setOffset(obsoffset_opt[0]); - imgReaderObs.setScale(obsscale_opt[0]); + if(obsoffset_opt.size()) + imgReaderObs.setOffset(obsoffset_opt[0]); + if(obsscale_opt.size()) + imgReaderObs.setScale(obsscale_opt[0]); //calculate regression between model and observation if(verbose_opt[0]) cout << "Calculating regression for " << imgReaderModel2.getFileName() << " " << imgReaderObs.getFileName() << endl; if(regSensor_opt[0]) - errObs=imgreg.getRMSE(imgReaderModel2,imgReaderObs,c0obs,c1obs,verbose_opt[0]); + errObs=imgreg.getRMSE(imgReaderModel2,imgReaderObs,c0obs,c1obs,0,0,verbose_opt[0]); else{ c0obs=0; c1obs=1; @@ -531,9 +641,11 @@ int main(int argc,char **argv) { } ImgReaderGdal imgReaderEst(input); imgReaderEst.setNoData(obsnodata_opt); - imgReaderEst.setOffset(obsoffset_opt[0]); - imgReaderEst.setScale(obsscale_opt[0]); - + if(obsoffset_opt.size()) + imgReaderEst.setOffset(obsoffset_opt[0]); + if(obsscale_opt.size()) + imgReaderEst.setScale(obsscale_opt[0]); + vector< vector<double> > obsLineVector(down_opt[0]); vector<double> obsLineBuffer; vector<double> obsWindowBuffer;//buffer for observation to calculate average corresponding to model pixel @@ -547,6 +659,7 @@ int main(int argc,char **argv) { vector<double> uncertReadBuffer; vector<double> estWriteBuffer(ncol); vector<double> uncertWriteBuffer(ncol); + // vector<double> lineMask; //initialize obsLineVector assert(down_opt[0]%2);//window size must be odd @@ -562,12 +675,12 @@ int main(int argc,char **argv) { imgReaderEst.readData(estReadBuffer,GDT_Float64,irow,0); imgReaderEst.readData(uncertReadBuffer,GDT_Float64,irow,1); //read model2 in case current estimate is nodata - imgReaderEst.image2geo(0,irow,x,y); - imgReaderModel2.geo2image(x,y,modCol,modRow); + imgReaderEst.image2geo(0,irow,geox,geoy); + imgReaderModel2.geo2image(geox,geoy,modCol,modRow); assert(modRow>=0&&modRow<imgReaderModel2.nrOfRow()); imgReaderModel2.readData(model2LineBuffer,GDT_Float64,modRow,0); - imgReaderModel1.geo2image(x,y,modCol,modRow); + imgReaderModel1.geo2image(geox,geoy,modCol,modRow); assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow()); imgReaderModel1.readData(model1LineBuffer,GDT_Float64,modRow,0); @@ -581,7 +694,37 @@ int main(int argc,char **argv) { if(imgReaderObs.nrOfBand()>1) imgReaderObs.readData(uncertObsLineBuffer,GDT_Float64,irow,1); } + // double oldRowMask=-1;//keep track of row mask to optimize number of line readings for(int icol=0;icol<imgWriterEst.nrOfCol();++icol){ + imgReaderEst.image2geo(icol,irow,geox,geoy); + // bool masked=false; + // if(mask_opt.size()){ + // //read mask + // maskReader.geo2image(geox,geoy,colMask,rowMask); + // colMask=static_cast<int>(colMask); + // rowMask=static_cast<int>(rowMask); + // if(rowMask>=0&&rowMask<maskReader.nrOfRow()&&colMask>=0&&colMask<maskReader.nrOfCol()){ + // if(static_cast<int>(rowMask)!=static_cast<int>(oldRowMask)){ + // try{ + // maskReader.readData(lineMask,GDT_Int16,static_cast<int>(rowMask)); + // } + // catch(string errorstring){ + // cerr << errorstring << endl; + // exit(1); + // } + // catch(...){ + // cerr << "error catched" << std::endl; + // exit(3); + // } + // oldRowMask=rowMask; + // } + // masked=(maskReader.isNoData(lineMask[colMask])); + // } + // estWriteBuffer[icol]=msknodata_opt[0]; + // uncertWriteBuffer[icol]=msknodata_opt[0]; + // continue;//next column + // } + int minCol=(icol>down_opt[0]/2) ? icol-down_opt[0]/2 : 0; int maxCol=(icol+down_opt[0]/2<imgReaderEst.nrOfCol()) ? icol+down_opt[0]/2 : imgReaderEst.nrOfCol()-1; int minRow=(irow>down_opt[0]/2) ? irow-down_opt[0]/2 : 0; @@ -600,8 +743,7 @@ int main(int argc,char **argv) { double estMeanValue=0;//stat.mean(estWindowBuffer); double nvalid=0; //time update - imgReaderEst.image2geo(icol,irow,x,y); - imgReaderModel2.geo2image(x,y,modCol,modRow); + imgReaderModel2.geo2image(geox,geoy,modCol,modRow); assert(modRow>=0&&modRow<imgReaderModel2.nrOfRow()); double modValue=model2LineBuffer[modCol]; bool estNodata=imgReaderEst.isNoData(estValue);//validity of current estimate @@ -619,21 +761,21 @@ int main(int argc,char **argv) { else{ if(window_opt[0]>0){ try{ - // imgReaderModel2.geo2image(x,y,modCol,modRow);//did that already + // imgReaderModel2.geo2image(geox,geoy,modCol,modRow);//did that already minCol=(modCol>window_opt[0]/2) ? modCol-window_opt[0]/2 : 0; maxCol=(modCol+window_opt[0]/2<imgReaderModel2.nrOfCol()) ? modCol+window_opt[0]/2 : imgReaderModel2.nrOfCol()-1; minRow=(modRow>window_opt[0]/2) ? modRow-window_opt[0]/2 : 0; maxRow=(modRow+window_opt[0]/2<imgReaderModel2.nrOfRow()) ? modRow+window_opt[0]/2 : imgReaderModel2.nrOfRow()-1; imgReaderModel2.readDataBlock(model2buffer,GDT_Float64,minCol,maxCol,minRow,maxRow,0); - imgReaderModel1.geo2image(x,y,modCol,modRow); + imgReaderModel1.geo2image(geox,geoy,modCol,modRow); assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow()); minCol=(modCol>window_opt[0]/2) ? modCol-window_opt[0]/2 : 0; maxCol=(modCol+window_opt[0]/2<imgReaderModel1.nrOfCol()) ? modCol+window_opt[0]/2 : imgReaderModel1.nrOfCol()-1; minRow=(modRow>window_opt[0]/2) ? modRow-window_opt[0]/2 : 0; maxRow=(modRow+window_opt[0]/2<imgReaderModel1.nrOfRow()) ? modRow+window_opt[0]/2 : imgReaderModel1.nrOfRow()-1; imgReaderModel1.readDataBlock(model1buffer,GDT_Float64,minCol,maxCol,minRow,maxRow,0); - // imgReaderEst.image2geo(icol,irow,x,y); + // imgReaderEst.image2geo(icol,irow,geox,geoy); } catch(string errorString){ cerr << "Error reading data block for " << minCol << "-" << maxCol << ", " << minRow << "-" << maxRow << endl; @@ -690,31 +832,29 @@ int main(int argc,char **argv) { if(update&&!imgReaderObs.isNoData(obsLineBuffer[icol])){ double kalmanGain=1; double uncertObs=uncertObs_opt[0]; - bool doUpdate=true; - if(deltaObs_opt.size()){ + if(uncertObsLineBuffer.size()>icol) + uncertObs=uncertObsLineBuffer[icol]; + else if(weight_opt.size()||deltaObs_opt.size()){ statfactory::StatFactory statobs; statobs.setNoDataValues(obsnodata_opt); double obsMeanValue=statobs.mean(obsWindowBuffer); double difference=(obsMeanValue-c0obs)/c1obs-modValue; - if(modValue){ - difference/=modValue;//relative difference - difference*=100;//in percentage - doUpdate=(difference>=deltaObs_opt[0]);//lower bound - doUpdate&=(difference<=deltaObs_opt[1]);//upper bound - //todo: use deltaObs to set observation uncertainty instead of setting doUpdate - if(-difference>uncertObs_opt[0]) - uncertObs=-difference; - } + if(modValue&&deltaObs_opt.size()){ + double relativeDifference=100.0*difference/modValue; + if(relativeDifference<deltaObs_opt[0])//lower bound + kalmanGain=0; + else if(relativeDifference>deltaObs_opt[1])//upper bound + kalmanGain=0; + } + uncertObs=-weight_opt[0]*difference; + if(uncertObs<=0) + uncertObs=0; } - if(doUpdate){ - if(uncertObsLineBuffer.size()>icol) - uncertObs=uncertObsLineBuffer[icol]; + if(kalmanGain>0){ if((uncertWriteBuffer[icol]+uncertObs)>eps_opt[0]) - kalmanGain=uncertWriteBuffer[icol]/(uncertWriteBuffer[icol]+uncertObs); - assert(kalmanGain<=1); + kalmanGain=uncertWriteBuffer[icol]/(uncertWriteBuffer[icol]+uncertObs); } - else - kalmanGain=0; + assert(kalmanGain<=1); estWriteBuffer[icol]+=kalmanGain*(obsLineBuffer[icol]-estWriteBuffer[icol]); uncertWriteBuffer[icol]*=(1-kalmanGain); } @@ -767,8 +907,10 @@ int main(int argc,char **argv) { try{ imgReaderModel1.open(model_opt.back()); imgReaderModel1.setNoData(modnodata_opt); - imgReaderModel1.setOffset(modoffset_opt[0]); - imgReaderModel1.setScale(modscale_opt[0]); + if(modoffset_opt.size()) + imgReaderModel1.setOffset(modoffset_opt[0]); + if(modscale_opt.size()) + imgReaderModel1.setScale(modscale_opt[0]); } catch(string errorString){ cerr << errorString << endl; @@ -783,8 +925,6 @@ int main(int argc,char **argv) { double minValue, maxValue, meanValue, stdDev; void* pProgressData; rasterBand->ComputeStatistics(0,&minValue,&maxValue,&meanValue,&stdDev,pfnProgress,pProgressData); - double x=0; - double y=0; double modRow=0; double modCol=0; if(relobsindex.back()<model_opt.size()-1){//initialize output_opt.back() as model[0] @@ -795,17 +935,46 @@ int main(int argc,char **argv) { vector<double> estReadBuffer; vector<double> estWriteBuffer(ncol); vector<double> uncertWriteBuffer(ncol); - imgWriterEst.image2geo(0,irow,x,y); - imgReaderModel1.geo2image(x,y,modCol,modRow); + // vector<double> lineMask; + imgWriterEst.image2geo(0,irow,geox,geoy); + imgReaderModel1.geo2image(geox,geoy,modCol,modRow); assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow()); try{ imgReaderModel1.readData(estReadBuffer,GDT_Float64,modRow); //simple nearest neighbor //stat.nearUp(estReadBuffer,estWriteBuffer); + // double oldRowMask=-1;//keep track of row mask to optimize number of line readings for(int icol=0;icol<imgWriterEst.nrOfCol();++icol){ - imgWriterEst.image2geo(icol,irow,x,y); - imgReaderModel1.geo2image(x,y,modCol,modRow); + imgWriterEst.image2geo(icol,irow,geox,geoy); + // bool masked=false; + // if(mask_opt.size()){ + // //read mask + // maskReader.geo2image(geox,geoy,colMask,rowMask); + // colMask=static_cast<int>(colMask); + // rowMask=static_cast<int>(rowMask); + // if(rowMask>=0&&rowMask<maskReader.nrOfRow()&&colMask>=0&&colMask<maskReader.nrOfCol()){ + // if(static_cast<int>(rowMask)!=static_cast<int>(oldRowMask)){ + // try{ + // maskReader.readData(lineMask,GDT_Int16,static_cast<int>(rowMask)); + // } + // catch(string errorstring){ + // cerr << errorstring << endl; + // exit(1); + // } + // catch(...){ + // cerr << "error catched" << std::endl; + // exit(3); + // } + // oldRowMask=rowMask; + // } + // masked=(maskReader.isNoData(lineMask[colMask])); + // } + // estWriteBuffer[icol]=msknodata_opt[0]; + // uncertWriteBuffer[icol]=msknodata_opt[0]; + // continue;//next column + // } + imgReaderModel1.geo2image(geox,geoy,modCol,modRow); estWriteBuffer[icol]=estReadBuffer[modCol]; uncertWriteBuffer[icol]=uncertModel_opt[0]*stdDev; } @@ -826,35 +995,68 @@ int main(int argc,char **argv) { imgReaderObs.open(observation_opt.back()); imgReaderObs.getGeoTransform(geotransform); imgReaderObs.setNoData(obsnodata_opt); - imgReaderObs.setOffset(obsoffset_opt[0]); - imgReaderObs.setScale(obsscale_opt[0]); - + if(obsoffset_opt.size()) + imgReaderObs.setOffset(obsoffset_opt[0]); + if(obsscale_opt.size()) + imgReaderObs.setScale(obsscale_opt[0]); + if(regSensor_opt[0]) - errObs=imgreg.getRMSE(imgReaderModel1,imgReaderObs,c0obs,c1obs,verbose_opt[0]); + errObs=imgreg.getRMSE(imgReaderModel1,imgReaderObs,c0obs,c1obs,0,0,verbose_opt[0]); else{ c0obs=0; c1obs=1; errObs=0; } + if(verbose_opt[0]) + cout << "c0obs, c1obs: " << c0obs << ", " << c1obs << endl; for(int irow=0;irow<nrow;++irow){ vector<double> estReadBuffer; - imgWriterEst.image2geo(0,irow,x,y); - imgReaderModel1.geo2image(x,y,modCol,modRow); + imgWriterEst.image2geo(0,irow,geox,geoy); + imgReaderModel1.geo2image(geox,geoy,modCol,modRow); assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow()); imgReaderModel1.readData(estReadBuffer,GDT_Float64,modRow); vector<double> obsLineBuffer; vector<double> estWriteBuffer(ncol); vector<double> uncertWriteBuffer(ncol); vector<double> uncertObsLineBuffer; + // vector<double> lineMask; // imgReaderObs.readData(estWriteBuffer,GDT_Float64,irow,0); imgReaderObs.readData(obsLineBuffer,GDT_Float64,irow,0); if(imgReaderObs.nrOfBand()>1) imgReaderObs.readData(uncertObsLineBuffer,GDT_Float64,irow,1); + // double oldRowMask=-1;//keep track of row mask to optimize number of line readings for(int icol=0;icol<imgWriterEst.nrOfCol();++icol){ - imgWriterEst.image2geo(icol,irow,x,y); - imgReaderModel1.geo2image(x,y,modCol,modRow); + imgWriterEst.image2geo(icol,irow,geox,geoy); + // bool masked=false; + // if(mask_opt.size()){ + // //read mask + // maskReader.geo2image(geox,geoy,colMask,rowMask); + // colMask=static_cast<int>(colMask); + // rowMask=static_cast<int>(rowMask); + // if(rowMask>=0&&rowMask<maskReader.nrOfRow()&&colMask>=0&&colMask<maskReader.nrOfCol()){ + // if(static_cast<int>(rowMask)!=static_cast<int>(oldRowMask)){ + // try{ + // maskReader.readData(lineMask,GDT_Int16,static_cast<int>(rowMask)); + // } + // catch(string errorstring){ + // cerr << errorstring << endl; + // exit(1); + // } + // catch(...){ + // cerr << "error catched" << std::endl; + // exit(3); + // } + // oldRowMask=rowMask; + // } + // masked=(maskReader.isNoData(lineMask[colMask])); + // } + // estWriteBuffer[icol]=msknodata_opt[0]; + // uncertWriteBuffer[icol]=msknodata_opt[0]; + // continue;//next column + // } + imgReaderModel1.geo2image(geox,geoy,modCol,modRow); assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow()); double modValue=estReadBuffer[modCol]; if(imgReaderModel1.isNoData(modValue)){//model is nodata: retain observation @@ -893,9 +1095,10 @@ int main(int argc,char **argv) { if(!imgReaderObs.isNoData(obsLineBuffer[icol])){ double kalmanGain=1; - bool doUpdate=true; double uncertObs=uncertObs_opt[0]; - if(deltaObs_opt.size()){ + if(uncertObsLineBuffer.size()>icol) + uncertObs=uncertObsLineBuffer[icol]; + else if(weight_opt.size()||deltaObs_opt.size()){ vector<double> obsWindowBuffer;//buffer for observation to calculate average corresponding to model pixel int minCol=(icol>down_opt[0]/2) ? icol-down_opt[0]/2 : 0; int maxCol=(icol+down_opt[0]/2<imgReaderObs.nrOfCol()) ? icol+down_opt[0]/2 : imgReaderObs.nrOfCol()-1; @@ -907,29 +1110,26 @@ int main(int argc,char **argv) { statobs.setNoDataValues(obsnodata_opt); double obsMeanValue=(statobs.mean(obsWindowBuffer)-c0obs)/c1obs; double difference=obsMeanValue-modValue; - if(modValue){ - difference/=modValue;//relative difference - difference*=100;//in percentage - doUpdate=(difference>=deltaObs_opt[0]);//lower bound - doUpdate&=(difference<=deltaObs_opt[1]);//upper bound - //todo: use deltaObs to set observation uncertainty instead of setting doUpdate - if(-difference>uncertObs_opt[0]) - uncertObs=-difference; - } - if(verbose_opt[0]>1){ - if(!doUpdate) - cout << "obsMeanValue:" << obsMeanValue << ", modValue: " << modValue << endl; - } + if(modValue&&deltaObs_opt.size()){ + double relativeDifference=100.0*difference/modValue; + if(relativeDifference<deltaObs_opt[0])//lower bound + kalmanGain=0; + else if(relativeDifference>deltaObs_opt[1])//upper bound + kalmanGain=0; + } + uncertObs=-weight_opt[0]*difference; + if(uncertObs<=0) + uncertObs=0; + if(verbose_opt[0]>1) + cout << "obsMeanValue:" << obsMeanValue << ", modValue: " << modValue << endl; } - if(doUpdate){ - if(uncertObsLineBuffer.size()>icol) - uncertObs=uncertObsLineBuffer[icol]; + if(kalmanGain>0){ if((uncertWriteBuffer[icol]+uncertObs)>eps_opt[0]) kalmanGain=uncertWriteBuffer[icol]/(uncertWriteBuffer[icol]+uncertObs); - assert(kalmanGain<=1); - estWriteBuffer[icol]+=kalmanGain*(obsLineBuffer[icol]-estWriteBuffer[icol]); - uncertWriteBuffer[icol]*=(1-kalmanGain); } + assert(kalmanGain<=1); + estWriteBuffer[icol]+=kalmanGain*(obsLineBuffer[icol]-estWriteBuffer[icol]); + uncertWriteBuffer[icol]*=(1-kalmanGain); } } imgWriterEst.writeData(estWriteBuffer,GDT_Float64,irow,0); @@ -968,12 +1168,16 @@ int main(int argc,char **argv) { //calculate regression between two subsequence model inputs imgReaderModel1.open(model_opt[modindex+1]); imgReaderModel1.setNoData(modnodata_opt); - imgReaderModel1.setOffset(modoffset_opt[0]); - imgReaderModel1.setScale(modscale_opt[0]); + if(modoffset_opt.size()) + imgReaderModel1.setOffset(modoffset_opt[0]); + if(modscale_opt.size()) + imgReaderModel1.setScale(modscale_opt[0]); imgReaderModel2.open(model_opt[modindex]); imgReaderModel2.setNoData(modnodata_opt); - imgReaderModel2.setOffset(modoffset_opt[0]); - imgReaderModel2.setScale(modscale_opt[0]); + if(modoffset_opt.size()) + imgReaderModel2.setOffset(modoffset_opt[0]); + if(modscale_opt.size()) + imgReaderModel2.setScale(modscale_opt[0]); //calculate regression //we could re-use the points from second image from last run, but //to keep it general, we must redo it (overlap might have changed) @@ -983,8 +1187,10 @@ int main(int argc,char **argv) { if(verbose_opt[0]) cout << "Calculating regression for " << imgReaderModel1.getFileName() << " " << imgReaderModel2.getFileName() << endl; - double errMod=imgreg.getRMSE(imgReaderModel1,imgReaderModel2,c0modGlobal,c1modGlobal); + double errMod=imgreg.getRMSE(imgReaderModel1,imgReaderModel2,c0modGlobal,c1modGlobal,0,0); // double errMod=imgreg.getRMSE(imgReaderModel1,imgReaderModel2,c0modGlobal,c1modGlobal,verbose_opt[0]); + if(verbose_opt[0]) + cout << "c0modGlobal, c1modGlobal: " << c0modGlobal << ", " << c1modGlobal << endl; bool update=false; if(obsindex<relobsindex.size()){ @@ -997,13 +1203,15 @@ int main(int argc,char **argv) { imgReaderObs.open(observation_opt[obsindex]); imgReaderObs.getGeoTransform(geotransform); imgReaderObs.setNoData(obsnodata_opt); - imgReaderObs.setOffset(obsoffset_opt[0]); - imgReaderObs.setScale(obsscale_opt[0]); + if(obsoffset_opt.size()) + imgReaderObs.setOffset(obsoffset_opt[0]); + if(obsscale_opt.size()) + imgReaderObs.setScale(obsscale_opt[0]); //calculate regression between model and observation if(verbose_opt[0]) cout << "Calculating regression for " << imgReaderModel2.getFileName() << " " << imgReaderObs.getFileName() << endl; if(regSensor_opt[0]) - errObs=imgreg.getRMSE(imgReaderModel2,imgReaderObs,c0obs,c1obs,verbose_opt[0]); + errObs=imgreg.getRMSE(imgReaderModel2,imgReaderObs,c0obs,c1obs,0,0,verbose_opt[0]); else{ c0obs=0; c1obs=1; @@ -1023,9 +1231,11 @@ int main(int argc,char **argv) { } ImgReaderGdal imgReaderEst(input); imgReaderEst.setNoData(obsnodata_opt); - imgReaderEst.setOffset(obsoffset_opt[0]); - imgReaderEst.setScale(obsscale_opt[0]); - + if(obsoffset_opt.size()) + imgReaderEst.setOffset(obsoffset_opt[0]); + if(obsscale_opt.size()) + imgReaderEst.setScale(obsscale_opt[0]); + vector< vector<double> > obsLineVector(down_opt[0]); vector<double> obsLineBuffer; vector<double> obsWindowBuffer;//buffer for observation to calculate average corresponding to model pixel @@ -1039,6 +1249,7 @@ int main(int argc,char **argv) { vector<double> uncertReadBuffer; vector<double> estWriteBuffer(ncol); vector<double> uncertWriteBuffer(ncol); + // vector<double> lineMask; //initialize obsLineVector assert(down_opt[0]%2);//window size must be odd @@ -1054,12 +1265,12 @@ int main(int argc,char **argv) { imgReaderEst.readData(estReadBuffer,GDT_Float64,irow,0); imgReaderEst.readData(uncertReadBuffer,GDT_Float64,irow,1); //read model2 in case current estimate is nodata - imgReaderEst.image2geo(0,irow,x,y); - imgReaderModel2.geo2image(x,y,modCol,modRow); + imgReaderEst.image2geo(0,irow,geox,geoy); + imgReaderModel2.geo2image(geox,geoy,modCol,modRow); assert(modRow>=0&&modRow<imgReaderModel2.nrOfRow()); imgReaderModel2.readData(model2LineBuffer,GDT_Float64,modRow,0); - imgReaderModel1.geo2image(x,y,modCol,modRow); + imgReaderModel1.geo2image(geox,geoy,modCol,modRow); assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow()); imgReaderModel1.readData(model1LineBuffer,GDT_Float64,modRow,0); @@ -1073,7 +1284,36 @@ int main(int argc,char **argv) { if(imgReaderObs.nrOfBand()>1) imgReaderObs.readData(uncertObsLineBuffer,GDT_Float64,irow,1); } + // double oldRowMask=-1;//keep track of row mask to optimize number of line readings for(int icol=0;icol<imgWriterEst.nrOfCol();++icol){ + imgReaderEst.image2geo(icol,irow,geox,geoy); + // bool masked=false; + // if(mask_opt.size()){ + // //read mask + // maskReader.geo2image(geox,geoy,colMask,rowMask); + // colMask=static_cast<int>(colMask); + // rowMask=static_cast<int>(rowMask); + // if(rowMask>=0&&rowMask<maskReader.nrOfRow()&&colMask>=0&&colMask<maskReader.nrOfCol()){ + // if(static_cast<int>(rowMask)!=static_cast<int>(oldRowMask)){ + // try{ + // maskReader.readData(lineMask,GDT_Int16,static_cast<int>(rowMask)); + // } + // catch(string errorstring){ + // cerr << errorstring << endl; + // exit(1); + // } + // catch(...){ + // cerr << "error catched" << std::endl; + // exit(3); + // } + // oldRowMask=rowMask; + // } + // masked=(maskReader.isNoData(lineMask[colMask])); + // } + // estWriteBuffer[icol]=msknodata_opt[0]; + // uncertWriteBuffer[icol]=msknodata_opt[0]; + // continue;//next column + // } int minCol=(icol>down_opt[0]/2) ? icol-down_opt[0]/2 : 0; int maxCol=(icol+down_opt[0]/2<imgReaderEst.nrOfCol()) ? icol+down_opt[0]/2 : imgReaderEst.nrOfCol()-1; int minRow=(irow>down_opt[0]/2) ? irow-down_opt[0]/2 : 0; @@ -1092,8 +1332,7 @@ int main(int argc,char **argv) { double estMeanValue=0;//stat.mean(estWindowBuffer); double nvalid=0; //time update - imgReaderEst.image2geo(icol,irow,x,y); - imgReaderModel2.geo2image(x,y,modCol,modRow); + imgReaderModel2.geo2image(geox,geoy,modCol,modRow); assert(modRow>=0&&modRow<imgReaderModel2.nrOfRow()); double modValue=model2LineBuffer[modCol]; bool estNodata=imgReaderEst.isNoData(estValue); @@ -1111,21 +1350,21 @@ int main(int argc,char **argv) { else{ if(window_opt[0]>0){ try{ - // imgReaderModel2.geo2image(x,y,modCol,modRow);//did that already + // imgReaderModel2.geo2image(geox,geoy,modCol,modRow);//did that already minCol=(modCol>window_opt[0]/2) ? modCol-window_opt[0]/2 : 0; maxCol=(modCol+window_opt[0]/2<imgReaderModel2.nrOfCol()) ? modCol+window_opt[0]/2 : imgReaderModel2.nrOfCol()-1; minRow=(modRow>window_opt[0]/2) ? modRow-window_opt[0]/2 : 0; maxRow=(modRow+window_opt[0]/2<imgReaderModel2.nrOfRow()) ? modRow+window_opt[0]/2 : imgReaderModel2.nrOfRow()-1; imgReaderModel2.readDataBlock(model2buffer,GDT_Float64,minCol,maxCol,minRow,maxRow,0); - imgReaderModel1.geo2image(x,y,modCol,modRow); + imgReaderModel1.geo2image(geox,geoy,modCol,modRow); assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow()); minCol=(modCol>window_opt[0]/2) ? modCol-window_opt[0]/2 : 0; maxCol=(modCol+window_opt[0]/2<imgReaderModel1.nrOfCol()) ? modCol+window_opt[0]/2 : imgReaderModel1.nrOfCol()-1; minRow=(modRow>window_opt[0]/2) ? modRow-window_opt[0]/2 : 0; maxRow=(modRow+window_opt[0]/2<imgReaderModel1.nrOfRow()) ? modRow+window_opt[0]/2 : imgReaderModel1.nrOfRow()-1; imgReaderModel1.readDataBlock(model1buffer,GDT_Float64,minCol,maxCol,minRow,maxRow,0); - // imgReaderEst.image2geo(icol,irow,x,y); + // imgReaderEst.image2geo(icol,irow,geox,geoy); } catch(string errorString){ cerr << "Error reading data block for " << minCol << "-" << maxCol << ", " << minRow << "-" << maxRow << endl; @@ -1182,31 +1421,29 @@ int main(int argc,char **argv) { if(update&&!imgReaderObs.isNoData(obsLineBuffer[icol])){ double kalmanGain=1; double uncertObs=uncertObs_opt[0]; - bool doUpdate=true; - if(deltaObs_opt.size()){ + if(uncertObsLineBuffer.size()>icol) + uncertObs=uncertObsLineBuffer[icol]; + else if(weight_opt.size()||deltaObs_opt.size()){ statfactory::StatFactory statobs; statobs.setNoDataValues(obsnodata_opt); double obsMeanValue=statobs.mean(obsWindowBuffer); double difference=(obsMeanValue-c0obs)/c1obs-modValue; - if(modValue){ - difference/=modValue;//relative difference - difference*=100;//in percentage - doUpdate=(difference>=deltaObs_opt[0]);//lower bound - doUpdate&=(difference<=deltaObs_opt[1]);//upper bound - //todo: use deltaObs to set observation uncertainty instead of setting doUpdate - if(-difference>uncertObs_opt[0]) - uncertObs=-difference; - } + if(modValue&&deltaObs_opt.size()){ + double relativeDifference=100.0*difference/modValue; + if(relativeDifference<deltaObs_opt[0])//lower bound + kalmanGain=0; + else if(relativeDifference>deltaObs_opt[1])//upper bound + kalmanGain=0; + } + uncertObs=-weight_opt[0]*difference; + if(uncertObs<=0) + uncertObs=0; } - if(doUpdate){ - if(uncertObsLineBuffer.size()>icol) - uncertObs=uncertObsLineBuffer[icol]; + if(kalmanGain>0){ if((uncertWriteBuffer[icol]+uncertObs)>eps_opt[0]) kalmanGain=uncertWriteBuffer[icol]/(uncertWriteBuffer[icol]+uncertObs); - assert(kalmanGain<=1); } - else - kalmanGain=0; + assert(kalmanGain<=1); estWriteBuffer[icol]+=kalmanGain*(obsLineBuffer[icol]-estWriteBuffer[icol]); uncertWriteBuffer[icol]*=(1-kalmanGain); } @@ -1277,12 +1514,16 @@ int main(int argc,char **argv) { ImgReaderGdal imgReaderForward(inputfw); ImgReaderGdal imgReaderBackward(inputbw); imgReaderForward.setNoData(obsnodata_opt); - imgReaderForward.setOffset(obsoffset_opt[0]); - imgReaderForward.setScale(obsscale_opt[0]); + if(obsoffset_opt.size()) + imgReaderForward.setOffset(obsoffset_opt[0]); + if(obsscale_opt.size()) + imgReaderForward.setScale(obsscale_opt[0]); imgReaderBackward.setNoData(obsnodata_opt); - imgReaderBackward.setOffset(obsoffset_opt[0]); - imgReaderBackward.setScale(obsscale_opt[0]); - + if(obsoffset_opt.size()) + imgReaderBackward.setOffset(obsoffset_opt[0]); + if(obsscale_opt.size()) + imgReaderBackward.setScale(obsscale_opt[0]); + vector<double> estForwardBuffer; vector<double> estBackwardBuffer; vector<double> uncertObsLineBuffer; @@ -1291,6 +1532,7 @@ int main(int argc,char **argv) { vector<double> uncertReadBuffer; vector<double> estWriteBuffer(ncol); vector<double> uncertWriteBuffer(ncol); + // vector<double> lineMask; bool update=false; if(obsindex<relobsindex.size()){ @@ -1303,8 +1545,10 @@ int main(int argc,char **argv) { imgReaderObs.open(observation_opt[obsindex]); imgReaderObs.getGeoTransform(geotransform); imgReaderObs.setNoData(obsnodata_opt); - imgReaderObs.setOffset(obsoffset_opt[0]); - imgReaderObs.setScale(obsscale_opt[0]); + if(obsoffset_opt.size()) + imgReaderObs.setOffset(obsoffset_opt[0]); + if(obsscale_opt.size()) + imgReaderObs.setScale(obsscale_opt[0]); //calculate regression between model and observation } @@ -1324,7 +1568,36 @@ int main(int argc,char **argv) { imgReaderObs.readData(uncertObsLineBuffer,GDT_Float64,irow,1); } + // double oldRowMask=-1;//keep track of row mask to optimize number of line readings for(int icol=0;icol<imgWriterEst.nrOfCol();++icol){ + imgWriterEst.image2geo(icol,irow,geox,geoy); + // bool masked=false; + // if(mask_opt.size()){ + // //read mask + // maskReader.geo2image(geox,geoy,colMask,rowMask); + // colMask=static_cast<int>(colMask); + // rowMask=static_cast<int>(rowMask); + // if(rowMask>=0&&rowMask<maskReader.nrOfRow()&&colMask>=0&&colMask<maskReader.nrOfCol()){ + // if(static_cast<int>(rowMask)!=static_cast<int>(oldRowMask)){ + // try{ + // maskReader.readData(lineMask,GDT_Int16,static_cast<int>(rowMask)); + // } + // catch(string errorstring){ + // cerr << errorstring << endl; + // exit(1); + // } + // catch(...){ + // cerr << "error catched" << std::endl; + // exit(3); + // } + // oldRowMask=rowMask; + // } + // masked=(maskReader.isNoData(lineMask[colMask])); + // } + // estWriteBuffer[icol]=msknodata_opt[0]; + // uncertWriteBuffer[icol]=msknodata_opt[0]; + // continue;//next column + // } double A=estForwardBuffer[icol]; double B=estBackwardBuffer[icol]; double C=uncertForwardBuffer[icol]*uncertForwardBuffer[icol]; @@ -1389,5 +1662,7 @@ int main(int argc,char **argv) { } } } + // if(mask_opt.size()) + // maskReader.close(); } diff --git a/src/apps/pkstat.cc b/src/apps/pkstat.cc new file mode 100644 index 0000000..5caa40a --- /dev/null +++ b/src/apps/pkstat.cc @@ -0,0 +1,825 @@ +/********************************************************************** +pkstat.cc: program to calculate basic statistics from raster dataset +Copyright (C) 2008-2015 Pieter Kempeneers + +This file is part of pktools + +pktools is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +pktools is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with pktools. If not, see <http://www.gnu.org/licenses/>. +***********************************************************************/ +#include <iostream> +#include <fstream> +#include <math.h> +#include "base/Optionpk.h" +#include "algorithms/StatFactory.h" +#include "algorithms/ImgRegression.h" +using namespace std; + +int main(int argc, char *argv[]) +{ + Optionpk<string> input_opt("i","input","name of the input raster dataset"); + Optionpk<unsigned short> band_opt("b","band","band(s) on which to calculate statistics",0); + Optionpk<bool> filename_opt("f", "filename", "Shows image filename ", false); + Optionpk<bool> stat_opt("stats", "statistics", "Shows basic statistics (min,max, mean and stdDev of the raster datasets)", false); + Optionpk<double> ulx_opt("ulx", "ulx", "Upper left x value bounding box"); + Optionpk<double> uly_opt("uly", "uly", "Upper left y value bounding box"); + Optionpk<double> lrx_opt("lrx", "lrx", "Lower right x value bounding box"); + Optionpk<double> lry_opt("lry", "lry", "Lower right y value bounding box"); + Optionpk<double> nodata_opt("nodata","nodata","Set nodata value(s)"); + Optionpk<short> down_opt("down", "down", "Down sampling factor (for raster sample datasets only). Can be used to create grid points", 1); + Optionpk<unsigned int> random_opt("rnd", "rnd", "generate random numbers", 0); + + // Optionpk<bool> transpose_opt("t","transpose","transpose output",false); + // Optionpk<std::string> randdist_opt("dist", "dist", "distribution for generating random numbers, see http://www.gn/software/gsl/manual/gsl-ref_toc.html#TOC320 (only uniform and Gaussian supported yet)", "gaussian"); + // Optionpk<double> randa_opt("rnda", "rnda", "first parameter for random distribution (mean value in case of Gaussian)", 0); + // Optionpk<double> randb_opt("rndb", "rndb", "second parameter for random distribution (standard deviation in case of Gaussian)", 1); + Optionpk<bool> mean_opt("mean","mean","calculate median",false); + Optionpk<bool> median_opt("median","median","calculate median",false); + Optionpk<bool> var_opt("var","var","calculate variance",false); + Optionpk<bool> skewness_opt("skew","skewness","calculate skewness",false); + Optionpk<bool> kurtosis_opt("kurt","kurtosis","calculate kurtosis",false); + Optionpk<bool> stdev_opt("stdev","stdev","calculate standard deviation",false); + Optionpk<bool> sum_opt("sum","sum","calculate sum of column",false); + Optionpk<bool> minmax_opt("mm","minmax","calculate minimum and maximum value",false); + Optionpk<bool> min_opt("min","min","calculate minimum value",false); + Optionpk<bool> max_opt("max","max","calculate maximum value",false); + Optionpk<double> src_min_opt("src_min","src_min","start reading source from this minimum value"); + Optionpk<double> src_max_opt("src_max","src_max","stop reading source from this maximum value"); + Optionpk<bool> histogram_opt("hist","hist","calculate histogram",false); + Optionpk<bool> histogram2d_opt("hist2d","hist2d","calculate 2-dimensional histogram based on two images",false); + Optionpk<short> nbin_opt("nbin","nbin","number of bins to calculate histogram"); + Optionpk<bool> relative_opt("rel","relative","use percentiles for histogram to calculate histogram",false); + Optionpk<bool> kde_opt("kde","kde","Use Kernel density estimation when producing histogram. The standard deviation is estimated based on Silverman's rule of thumb",false); + Optionpk<bool> correlation_opt("cor","correlation","calculate Pearson produc-moment correlation coefficient between two raster datasets (defined by -c <col1> -c <col2>",false); + Optionpk<bool> rmse_opt("rmse","rmse","calculate root mean square error between two raster datasets",false); + Optionpk<bool> reg_opt("reg","regression","calculate linear regression between two raster datasets and get correlation coefficient",false); + Optionpk<bool> regerr_opt("regerr","regerr","calculate linear regression between two raster datasets and get root mean square error",false); + Optionpk<bool> preg_opt("preg","preg","calculate perpendicular regression between two raster datasets and get correlation coefficient",false); + Optionpk<bool> pregerr_opt("pregerr","pregerr","calculate perpendicular regression between two raster datasets and get root mean square error",false); + Optionpk<short> verbose_opt("v", "verbose", "verbose mode when positive", 0,2); + + src_min_opt.setHide(1); + src_max_opt.setHide(1); + // range_opt.setHide(1); + // transpose_opt.setHide(1); + + bool doProcess;//stop process when program was invoked with help option (-h --help) + try{ + //mandatory options + doProcess=input_opt.retrieveOption(argc,argv); + //optional options + band_opt.retrieveOption(argc,argv); + filename_opt.retrieveOption(argc,argv); + stat_opt.retrieveOption(argc,argv); + ulx_opt.retrieveOption(argc,argv); + uly_opt.retrieveOption(argc,argv); + lrx_opt.retrieveOption(argc,argv); + lry_opt.retrieveOption(argc,argv); + nodata_opt.retrieveOption(argc,argv); + down_opt.retrieveOption(argc,argv); + random_opt.retrieveOption(argc,argv); + // randdist_opt.retrieveOption(argc,argv); + // randa_opt.retrieveOption(argc,argv); + // randb_opt.retrieveOption(argc,argv); + mean_opt.retrieveOption(argc,argv); + median_opt.retrieveOption(argc,argv); + var_opt.retrieveOption(argc,argv); + stdev_opt.retrieveOption(argc,argv); + // skewness_opt.retrieveOption(argc,argv); + // kurtosis_opt.retrieveOption(argc,argv); + // sum_opt.retrieveOption(argc,argv); + minmax_opt.retrieveOption(argc,argv); + min_opt.retrieveOption(argc,argv); + max_opt.retrieveOption(argc,argv); + histogram_opt.retrieveOption(argc,argv); + nbin_opt.retrieveOption(argc,argv); + relative_opt.retrieveOption(argc,argv); + kde_opt.retrieveOption(argc,argv); + histogram2d_opt.retrieveOption(argc,argv); + correlation_opt.retrieveOption(argc,argv); + rmse_opt.retrieveOption(argc,argv); + reg_opt.retrieveOption(argc,argv); + regerr_opt.retrieveOption(argc,argv); + preg_opt.retrieveOption(argc,argv); + pregerr_opt.retrieveOption(argc,argv); + //advanced options + src_min_opt.retrieveOption(argc,argv); + src_max_opt.retrieveOption(argc,argv); + // range_opt.retrieveOption(argc,argv); + // transpose_opt.retrieveOption(argc,argv); + // comment_opt.retrieveOption(argc,argv); + verbose_opt.retrieveOption(argc,argv); + } + catch(string predefinedString){ + std::cout << predefinedString << std::endl; + exit(0); + } + if(!doProcess){ + cout << endl; + cout << "Usage: pkstat -i input [-c column]*" << endl; + cout << endl; + std::cout << "short option -h shows basic options only, use long option --help to show all options" << std::endl; + exit(0);//help was invoked, stop processing + } + + if(src_min_opt.size()){ + while(src_min_opt.size()<band_opt.size()) + src_min_opt.push_back(src_min_opt[0]); + } + if(src_max_opt.size()){ + while(src_max_opt.size()<band_opt.size()) + src_max_opt.push_back(src_max_opt[0]); + } + + unsigned int nbin=0; + double minX=0; + double minY=0; + double maxX=0; + double maxY=0; + double minValue=0; + double maxValue=0; + double meanValue=0; + double stdDev=0; + + const char* pszMessage; + void* pProgressArg=NULL; + GDALProgressFunc pfnProgress=GDALTermProgress; + double progress=0; + srand(time(NULL)); + + statfactory::StatFactory stat; + imgregression::ImgRegression imgreg; + + ImgReaderGdal imgReader; + if(input_opt.empty()){ + std::cerr << "No image dataset provided (use option -i). Use --help for help information"; + exit(0); + } + for(int ifile=0;ifile<input_opt.size();++ifile){ + try{ + imgReader.open(input_opt[ifile]); + } + catch(std::string errorstring){ + std::cout << errorstring << std::endl; + exit(0); + } + for(int inodata=0;inodata<nodata_opt.size();++inodata){ + if(!inodata) + imgReader.GDALSetNoDataValue(nodata_opt[0],band_opt[0]);//only single no data can be set in GDALRasterBand (used for ComputeStatistics) + imgReader.pushNoDataValue(nodata_opt[inodata]); + } + if(filename_opt[0]) + std::cout << " --input " << input_opt[ifile] << " "; + int nband=band_opt.size(); + for(int iband=0;iband<nband;++iband){ + if(stat_opt[0]||mean_opt[0]||var_opt[0]||stdev_opt[0]){ + assert(band_opt[iband]<imgReader.nrOfBand()); + GDALProgressFunc pfnProgress; + void* pProgressData; + GDALRasterBand* rasterBand; + + rasterBand=imgReader.getRasterBand(band_opt[iband]); + rasterBand->ComputeStatistics(0,&minValue,&maxValue,&meanValue,&stdDev,pfnProgress,pProgressData); + if(mean_opt[0]) + std::cout << "--mean " << meanValue << " "; + if(stdev_opt[0]) + std::cout << "--stdDev " << stdDev << " "; + if(var_opt[0]) + std::cout << "--var " << stdDev*stdDev << " "; + if(stat_opt[0]) + std::cout << "-min " << minValue << " -max " << maxValue << " --mean " << meanValue << " --stdDev " << stdDev << " "; + } + + if(minmax_opt[0]||min_opt[0]||max_opt[0]){ + assert(band_opt[iband]<imgReader.nrOfBand()); + if((ulx_opt.size()||uly_opt.size()||lrx_opt.size()||lry_opt.size())&&(imgReader.covers(ulx_opt[0],uly_opt[0],lrx_opt[0],lry_opt[0]))){ + double uli,ulj,lri,lrj; + imgReader.geo2image(ulx_opt[0],uly_opt[0],uli,ulj); + imgReader.geo2image(lrx_opt[0],lry_opt[0],lri,lrj); + imgReader.getMinMax(static_cast<int>(uli),static_cast<int>(lri),static_cast<int>(ulj),static_cast<int>(lrj),band_opt[iband],minValue,maxValue); + } + else + imgReader.getMinMax(minValue,maxValue,band_opt[iband],true); + if(minmax_opt[0]) + std::cout << "-min " << minValue << " -max " << maxValue << " "; + else{ + if(min_opt[0]) + std::cout << "-min " << minValue << " "; + if(max_opt[0]) + std::cout << "-max " << maxValue << " "; + } + } + } + if(histogram_opt[0]){//only for first selected band + assert(band_opt[0]<imgReader.nrOfBand()); + nbin=(nbin_opt.size())? nbin_opt[0]:0; + + imgReader.getMinMax(minValue,maxValue,band_opt[0]); + if(src_min_opt.size()) + minValue=src_min_opt[0]; + if(src_max_opt.size()) + maxValue=src_max_opt[0]; + if(minValue>=maxValue) + imgReader.getMinMax(minValue,maxValue,band_opt[0]); + + if(verbose_opt[0]) + cout << "number of valid pixels in image: " << imgReader.getNvalid(band_opt[0]) << endl; + std::vector<double> output; + double nsample=imgReader.getHistogram(output,minValue,maxValue,nbin,band_opt[0],kde_opt[0]); + + std::cout.precision(10); + for(int bin=0;bin<nbin;++bin){ + double binValue=0; + if(nbin==maxValue-minValue+1) + binValue=minValue+bin; + else + binValue=minValue+static_cast<double>(maxValue-minValue)*(bin+0.5)/nbin; + std::cout << binValue << " "; + if(relative_opt[0]||kde_opt[0]) + std::cout << 100.0*static_cast<double>(output[bin])/static_cast<double>(nsample) << std::endl; + else + std::cout << static_cast<double>(output[bin]) << std::endl; + } + } + if(histogram2d_opt[0]){ + if(band_opt.size()<2) + continue; + imgReader.getMinMax(minX,maxX,band_opt[0]); + imgReader.getMinMax(minY,maxY,band_opt[1]); + if(src_min_opt.size()){ + minX=src_min_opt[0]; + minY=src_min_opt[1]; + } + if(src_max_opt.size()){ + maxX=src_max_opt[0]; + maxY=src_max_opt[1]; + } + nbin=(nbin_opt.size())? nbin_opt[0]:0; + if(nbin<=1){ + std::cerr << "Warning: number of bins not defined, calculating bins from min and max value" << std::endl; + if(minX>=maxX) + imgReader.getMinMax(minX,maxX,band_opt[0]); + if(minY>=maxY) + imgReader.getMinMax(minY,maxY,band_opt[1]); + + minValue=(minX<minY)? minX:minY; + maxValue=(maxX>maxY)? maxX:maxY; + if(verbose_opt[0]) + std::cout << "min and max values: " << minValue << ", " << maxValue << std::endl; + nbin=maxValue-minValue+1; + } + assert(nbin>1); + double sigma=0; + //kernel density estimation as in http://en.wikipedia.org/wiki/Kernel_density_estimation + if(kde_opt[0]){ + assert(band_opt[0]<imgReader.nrOfBand()); + assert(band_opt[1]<imgReader.nrOfBand()); + GDALProgressFunc pfnProgress; + void* pProgressData; + GDALRasterBand* rasterBand; + double stdDev1=0; + double stdDev2=0; + rasterBand=imgReader.getRasterBand(band_opt[0]); + rasterBand->ComputeStatistics(0,&minValue,&maxValue,&meanValue,&stdDev1,pfnProgress,pProgressData); + rasterBand=imgReader.getRasterBand(band_opt[1]); + rasterBand->ComputeStatistics(0,&minValue,&maxValue,&meanValue,&stdDev2,pfnProgress,pProgressData); + + double estimatedSize=1.0*imgReader.getNvalid(band_opt[0])/down_opt[0]/down_opt[0]; + if(random_opt[0]>0) + estimatedSize*=random_opt[0]/100.0; + sigma=1.06*sqrt(stdDev1*stdDev2)*pow(estimatedSize,-0.2); + } + assert(nbin); + if(verbose_opt[0]){ + if(sigma>0) + std::cout << "calculating 2d kernel density estimate with sigma " << sigma << " for bands " << band_opt[0] << " and " << band_opt[1] << std::endl; + else + std::cout << "calculating 2d histogram for bands " << band_opt[0] << " and " << band_opt[1] << std::endl; + std::cout << "nbin: " << nbin << std::endl; + } + + + vector< vector<double> > output; + + if(maxX<=minX) + imgReader.getMinMax(minX,maxX,band_opt[0]); + if(maxY<=minY) + imgReader.getMinMax(minY,maxY,band_opt[1]); + + if(maxX<=minX){ + std::ostringstream s; + s<<"Error: could not calculate distribution (minX>=maxX)"; + throw(s.str()); + } + if(maxY<=minY){ + std::ostringstream s; + s<<"Error: could not calculate distribution (minY>=maxY)"; + throw(s.str()); + } + output.resize(nbin); + for(int i=0;i<nbin;++i){ + output[i].resize(nbin); + for(int j=0;j<nbin;++j) + output[i][j]=0; + } + int binX=0; + int binY=0; + vector<double> inputX(imgReader.nrOfCol()); + vector<double> inputY(imgReader.nrOfCol()); + unsigned long int nvalid=0; + for(int irow=0;irow<imgReader.nrOfRow();++irow){ + if(irow%down_opt[0]) + continue; + imgReader.readData(inputX,GDT_Float64,irow,band_opt[0]); + imgReader.readData(inputY,GDT_Float64,irow,band_opt[1]); + for(int icol=0;icol<imgReader.nrOfCol();++icol){ + if(icol%down_opt[0]) + continue; + if(random_opt[0]>0){ + double p=static_cast<double>(rand())/(RAND_MAX); + p*=100.0; + if(p>random_opt[0]) + continue;//do not select for now, go to next column + } + if(imgReader.isNoData(inputX[icol])) + continue; + if(imgReader.isNoData(inputY[icol])) + continue; + ++nvalid; + if(inputX[icol]>=maxX) + binX=nbin-1; + else if(inputX[icol]<=minX) + binX=0; + else + binX=static_cast<int>(static_cast<double>(inputX[icol]-minX)/(maxX-minX)*nbin); + if(inputY[icol]>=maxY) + binY=nbin-1; + else if(inputY[icol]<=minX) + binY=0; + else + binY=static_cast<int>(static_cast<double>(inputY[icol]-minY)/(maxY-minY)*nbin); + assert(binX>=0); + assert(binX<output.size()); + assert(binY>=0); + assert(binY<output[binX].size()); + if(sigma>0){ + //create kde for Gaussian basis function + //todo: speed up by calculating first and last bin with non-zero contriubtion... + for(int ibinX=0;ibinX<nbin;++ibinX){ + double centerX=minX+static_cast<double>(maxX-minX)*ibinX/nbin; + double pdfX=gsl_ran_gaussian_pdf(inputX[icol]-centerX, sigma); + for(int ibinY=0;ibinY<nbin;++ibinY){ + //calculate \integral_ibinX^(ibinX+1) + double centerY=minY+static_cast<double>(maxY-minY)*ibinY/nbin; + double pdfY=gsl_ran_gaussian_pdf(inputY[icol]-centerY, sigma); + output[ibinX][binY]+=pdfX*pdfY; + } + } + } + else + ++output[binX][binY]; + } + } + if(verbose_opt[0]) + cout << "number of valid pixels: " << nvalid << endl; + + for(int binX=0;binX<nbin;++binX){ + cout << endl; + for(int binY=0;binY<nbin;++binY){ + double binValueX=0; + if(nbin==maxX-minX+1) + binValueX=minX+binX; + else + binValueX=minX+static_cast<double>(maxX-minX)*(binX+0.5)/nbin; + double binValueY=0; + if(nbin==maxY-minY+1) + binValueY=minY+binY; + else + binValueY=minY+static_cast<double>(maxY-minY)*(binY+0.5)/nbin; + + double value=static_cast<double>(output[binX][binY]); + + if(relative_opt[0]) + value*=100.0/nvalid; + + cout << binValueX << " " << binValueY << " " << value << std::endl; + // double value=static_cast<double>(output[binX][binY])/nvalid; + // cout << (maxX-minX)*bin/(nbin-1)+minX << " " << (maxY-minY)*bin/(nbin-1)+minY << " " << value << std::endl; + } + } + } + if(reg_opt[0]){ + if(band_opt.size()<2) + continue; + imgreg.setDown(down_opt[0]); + imgreg.setThreshold(random_opt[0]); + double c0=0;//offset + double c1=1;//scale + double r2=imgreg.getR2(imgReader,band_opt[0],band_opt[1],c0,c1,verbose_opt[0]); + std::cout << "-c0 " << c0 << " -c1 " << c1 << " -r2 " << r2 << std::endl; + } + if(regerr_opt[0]){ + if(band_opt.size()<2) + continue; + imgreg.setDown(down_opt[0]); + imgreg.setThreshold(random_opt[0]); + double c0=0;//offset + double c1=1;//scale + double err=imgreg.getRMSE(imgReader,band_opt[0],band_opt[1],c0,c1,verbose_opt[0]); + std::cout << "-c0 " << c0 << " -c1 " << c1 << " -rmse " << err << std::endl; + } + if(rmse_opt[0]){ + if(band_opt.size()<2) + continue; + imgreg.setDown(down_opt[0]); + imgreg.setThreshold(random_opt[0]); + double c0=0;//offset + double c1=1;//scale + double err=imgreg.getRMSE(imgReader,band_opt[0],band_opt[1],c0,c1,verbose_opt[0]); + std::cout << " -rmse " << err << std::endl; + } + if(preg_opt[0]){ + if(band_opt.size()<2) + continue; + imgreg.setDown(down_opt[0]); + imgreg.setThreshold(random_opt[0]); + double c0=0;//offset + double c1=1;//scale + double r2=imgreg.pgetR2(imgReader,band_opt[0],band_opt[1],c0,c1,verbose_opt[0]); + std::cout << "-c0 " << c0 << " -c1 " << c1 << " -r2 " << r2 << std::endl; + } + imgReader.close(); + } + if(reg_opt[0]&&(input_opt.size()>1)){ + while(band_opt.size()<input_opt.size()) + band_opt.push_back(band_opt[0]); + imgreg.setDown(down_opt[0]); + imgreg.setThreshold(random_opt[0]); + double c0=0;//offset + double c1=1;//scale + ImgReaderGdal imgReader1(input_opt[0]); + ImgReaderGdal imgReader2(input_opt[1]); + double r2=imgreg.getR2(imgReader1,imgReader2,c0,c1,band_opt[0],band_opt[1],verbose_opt[0]); + std::cout << "-c0 " << c0 << " -c1 " << c1 << " -r2 " << r2 << std::endl; + imgReader1.close(); + imgReader2.close(); + } + if(preg_opt[0]&&(input_opt.size()>1)){ + while(band_opt.size()<input_opt.size()) + band_opt.push_back(band_opt[0]); + imgreg.setDown(down_opt[0]); + imgreg.setThreshold(random_opt[0]); + double c0=0;//offset + double c1=1;//scale + ImgReaderGdal imgReader1(input_opt[0]); + ImgReaderGdal imgReader2(input_opt[1]); + double r2=imgreg.pgetR2(imgReader1,imgReader2,c0,c1,band_opt[0],band_opt[1],verbose_opt[0]); + std::cout << "-c0 " << c0 << " -c1 " << c1 << " -r2 " << r2 << std::endl; + imgReader1.close(); + imgReader2.close(); + } + if(regerr_opt[0]&&(input_opt.size()>1)){ + while(band_opt.size()<input_opt.size()) + band_opt.push_back(band_opt[0]); + imgreg.setDown(down_opt[0]); + imgreg.setThreshold(random_opt[0]); + double c0=0;//offset + double c1=1;//scale + ImgReaderGdal imgReader1(input_opt[0]); + ImgReaderGdal imgReader2(input_opt[1]); + double err=imgreg.getRMSE(imgReader1,imgReader2,c0,c1,band_opt[0],band_opt[1],verbose_opt[0]); + std::cout << "-c0 " << c0 << " -c1 " << c1 << " -rmse " << err << std::endl; + imgReader1.close(); + imgReader2.close(); + } + if(rmse_opt[0]&&(input_opt.size()>1)){ + while(band_opt.size()<input_opt.size()) + band_opt.push_back(band_opt[0]); + imgreg.setDown(down_opt[0]); + imgreg.setThreshold(random_opt[0]); + double c0=0;//offset + double c1=1;//scale + ImgReaderGdal imgReader1(input_opt[0]); + ImgReaderGdal imgReader2(input_opt[1]); + double err=imgreg.getRMSE(imgReader1,imgReader2,c0,c1,band_opt[0],band_opt[1],verbose_opt[0]); + std::cout << "-rmse " << err << std::endl; + imgReader1.close(); + imgReader2.close(); + } + if(histogram2d_opt[0]&&(input_opt.size()>1)){ + imgReader.getMinMax(minX,maxX,band_opt[0]); + imgReader.getMinMax(minY,maxY,band_opt[1]); + if(src_min_opt.size()){ + while(src_min_opt.size()<input_opt.size()) + src_min_opt.push_back(src_min_opt[0]); + } + if(src_max_opt.size()){ + while(src_max_opt.size()<input_opt.size()) + src_max_opt.push_back(src_max_opt[0]); + } + if(src_min_opt.size()){ + minX=src_min_opt[0]; + minY=src_min_opt[1]; + } + if(src_max_opt.size()){ + maxX=src_max_opt[0]; + maxY=src_max_opt[1]; + } + nbin=(nbin_opt.size())? nbin_opt[0]:0; + ImgReaderGdal imgReader1(input_opt[0]); + ImgReaderGdal imgReader2(input_opt[1]); + for(int inodata=0;inodata<nodata_opt.size();++inodata){ + if(!inodata){ + imgReader1.GDALSetNoDataValue(nodata_opt[0],band_opt[0]);//only single no data can be set in GDALRasterBand (used for ComputeStatistics) + imgReader2.GDALSetNoDataValue(nodata_opt[0],band_opt[0]);//only single no data can be set in GDALRasterBand (used for ComputeStatistics) + } + imgReader1.pushNoDataValue(nodata_opt[inodata]); + imgReader2.pushNoDataValue(nodata_opt[inodata]); + } + if(nbin<=1){ + std::cerr << "Warning: number of bins not defined, calculating bins from min and max value" << std::endl; + // imgReader1.getMinMax(minX,maxX,band_opt[0]); + // imgReader2.getMinMax(minY,maxY,band_opt[0]); + if(minX>=maxX) + imgReader1.getMinMax(minX,maxX,band_opt[0]); + if(minY>=maxY) + imgReader2.getMinMax(minY,maxY,band_opt[1]); + + minValue=(minX<minY)? minX:minY; + maxValue=(maxX>maxY)? maxX:maxY; + if(verbose_opt[0]) + std::cout << "min and max values: " << minValue << ", " << maxValue << std::endl; + nbin=maxValue-minValue+1; + } + assert(nbin>1); + double sigma=0; + //kernel density estimation as in http://en.wikipedia.org/wiki/Kernel_density_estimation + if(kde_opt[0]){ + GDALProgressFunc pfnProgress; + void* pProgressData; + GDALRasterBand* rasterBand; + double stdDev1=0; + double stdDev2=0; + rasterBand=imgReader1.getRasterBand(band_opt[0]); + rasterBand->ComputeStatistics(0,&minValue,&maxValue,&meanValue,&stdDev1,pfnProgress,pProgressData); + rasterBand=imgReader2.getRasterBand(band_opt[0]); + rasterBand->ComputeStatistics(0,&minValue,&maxValue,&meanValue,&stdDev2,pfnProgress,pProgressData); + + //todo: think of smarter way how to estimate size (nodata!) + double estimatedSize=1.0*imgReader.getNvalid(band_opt[0])/down_opt[0]/down_opt[0]; + if(random_opt[0]>0) + estimatedSize*=random_opt[0]/100.0; + sigma=1.06*sqrt(stdDev1*stdDev2)*pow(estimatedSize,-0.2); + } + assert(nbin); + if(verbose_opt[0]){ + if(sigma>0) + std::cout << "calculating 2d kernel density estimate with sigma " << sigma << " for datasets " << input_opt[0] << " and " << input_opt[1] << std::endl; + else + std::cout << "calculating 2d histogram for datasets " << input_opt[0] << " and " << input_opt[1] << std::endl; + std::cout << "nbin: " << nbin << std::endl; + } + + vector< vector<double> > output; + + if(maxX<=minX) + imgReader1.getMinMax(minX,maxX,band_opt[0]); + if(maxY<=minY) + imgReader2.getMinMax(minY,maxY,band_opt[0]); + + if(maxX<=minX){ + std::ostringstream s; + s<<"Error: could not calculate distribution (minX>=maxX)"; + throw(s.str()); + } + if(maxY<=minY){ + std::ostringstream s; + s<<"Error: could not calculate distribution (minY>=maxY)"; + throw(s.str()); + } + if(verbose_opt[0]){ + cout << "minX: " << minX << endl; + cout << "maxX: " << maxX << endl; + cout << "minY: " << minY << endl; + cout << "maxY: " << maxY << endl; + } + output.resize(nbin); + for(int i=0;i<nbin;++i){ + output[i].resize(nbin); + for(int j=0;j<nbin;++j) + output[i][j]=0; + } + int binX=0; + int binY=0; + vector<double> inputX(imgReader1.nrOfCol()); + vector<double> inputY(imgReader2.nrOfCol()); + double nvalid=0; + double geoX=0; + double geoY=0; + double icol1=0; + double irow1=0; + double icol2=0; + double irow2=0; + for(int irow=0;irow<imgReader1.nrOfRow();++irow){ + if(irow%down_opt[0]) + continue; + irow1=irow; + imgReader1.image2geo(icol1,irow1,geoX,geoY); + imgReader2.geo2image(geoX,geoY,icol2,irow2); + irow2=static_cast<int>(irow2); + imgReader1.readData(inputX,GDT_Float64,irow1,band_opt[0]); + imgReader2.readData(inputY,GDT_Float64,irow2,band_opt[0]); + for(int icol=0;icol<imgReader.nrOfCol();++icol){ + if(icol%down_opt[0]) + continue; + icol1=icol; + if(random_opt[0]>0){ + double p=static_cast<double>(rand())/(RAND_MAX); + p*=100.0; + if(p>random_opt[0]) + continue;//do not select for now, go to next column + } + if(imgReader1.isNoData(inputX[icol])) + continue; + imgReader1.image2geo(icol1,irow1,geoX,geoY); + imgReader2.geo2image(geoX,geoY,icol2,irow2); + icol2=static_cast<int>(icol2); + if(imgReader2.isNoData(inputY[icol2])) + continue; + // ++nvalid; + if(inputX[icol1]>=maxX) + binX=nbin-1; + else if(inputX[icol]<=minX) + binX=0; + else + binX=static_cast<int>(static_cast<double>(inputX[icol1]-minX)/(maxX-minX)*nbin); + if(inputY[icol2]>=maxY) + binY=nbin-1; + else if(inputY[icol2]<=minY) + binY=0; + else + binY=static_cast<int>(static_cast<double>(inputY[icol2]-minY)/(maxY-minY)*nbin); + assert(binX>=0); + assert(binX<output.size()); + assert(binY>=0); + assert(binY<output[binX].size()); + if(sigma>0){ + //create kde for Gaussian basis function + //todo: speed up by calculating first and last bin with non-zero contriubtion... + for(int ibinX=0;ibinX<nbin;++ibinX){ + double centerX=minX+static_cast<double>(maxX-minX)*ibinX/nbin; + double pdfX=gsl_ran_gaussian_pdf(inputX[icol1]-centerX, sigma); + for(int ibinY=0;ibinY<nbin;++ibinY){ + //calculate \integral_ibinX^(ibinX+1) + double centerY=minY+static_cast<double>(maxY-minY)*ibinY/nbin; + double pdfY=gsl_ran_gaussian_pdf(inputY[icol2]-centerY, sigma); + output[ibinX][binY]+=pdfX*pdfY; + nvalid+=pdfX*pdfY; + } + } + } + else{ + ++output[binX][binY]; + ++nvalid; + } + } + } + if(verbose_opt[0]) + cout << "number of valid pixels: " << nvalid << endl; + for(int binX=0;binX<nbin;++binX){ + cout << endl; + for(int binY=0;binY<nbin;++binY){ + double binValueX=0; + if(nbin==maxX-minX+1) + binValueX=minX+binX; + else + binValueX=minX+static_cast<double>(maxX-minX)*(binX+0.5)/nbin; + double binValueY=0; + if(nbin==maxY-minY+1) + binValueY=minY+binY; + else + binValueY=minY+static_cast<double>(maxY-minY)*(binY+0.5)/nbin; + double value=static_cast<double>(output[binX][binY]); + + if(relative_opt[0]||kde_opt[0]) + value*=100.0/nvalid; + + cout << binValueX << " " << binValueY << " " << value << std::endl; + // double value=static_cast<double>(output[binX][binY])/nvalid; + // cout << (maxX-minX)*bin/(nbin-1)+minX << " " << (maxY-minY)*bin/(nbin-1)+minY << " " << value << std::endl; + } + } + imgReader1.close(); + imgReader2.close(); + } + + if(!histogram_opt[0]||histogram2d_opt[0]) + std::cout << std::endl; +} + +// int nband=(band_opt.size()) ? band_opt.size() : imgReader.nrOfBand(); + +// const char* pszMessage; +// void* pProgressArg=NULL; +// GDALProgressFunc pfnProgress=GDALTermProgress; +// double progress=0; +// srand(time(NULL)); + + +// statfactory::StatFactory stat; +// imgregression::ImgRegression imgreg; + +// pfnProgress(progress,pszMessage,pProgressArg); +// for(irow=0;irow<classReader.nrOfRow();++irow){ +// if(irow%down_opt[0]) +// continue; +// // classReader.readData(classBuffer,GDT_Int32,irow); +// classReader.readData(classBuffer,GDT_Float64,irow); +// double x,y;//geo coordinates +// double iimg,jimg;//image coordinates in img image +// for(icol=0;icol<classReader.nrOfCol();++icol){ +// if(icol%down_opt[0]) + // continue; + + + // if(rand_opt[0]>0){ + // gsl_rng* r=stat.getRandomGenerator(time(NULL)); + // //todo: init random number generator using time... + // if(verbose_opt[0]) + // std::cout << "generating " << rand_opt[0] << " random numbers: " << std::endl; + // for(unsigned int i=0;i<rand_opt[0];++i) + // std::cout << i << " " << stat.getRandomValue(r,randdist_opt[0],randa_opt[0],randb_opt[0]) << std::endl; + // } + + // imgreg.setDown(down_opt[0]); + // imgreg.setThreshold(threshold_opt[0]); + // double c0=0;//offset + // double c1=1;//scale + // double err=uncertNodata_opt[0];//start with high initial value in case we do not have first ob err=imgreg.getRMSE(imgReaderModel1,imgReader,c0,c1,verbose_opt[0]); + + // int nband=band_opt.size(); + // if(band_opt[0]<0) + // nband=imgReader.nrOfBand(); + // for(int iband=0;iband<nband;++iband){ + // unsigned short band_opt[iband]=(band_opt[0]<0)? iband : band_opt[iband]; + + // if(minmax_opt[0]||min_opt[0]||max_opt[0]){ + // assert(band_opt[iband]<imgReader.nrOfBand()); + // if((ulx_opt.size()||uly_opt.size()||lrx_opt.size()||lry_opt.size())&&(imgReader.covers(ulx_opt[0],uly_opt[0],lrx_opt[0],lry_opt[0]))){ + // double uli,ulj,lri,lrj; + // imgReader.geo2image(ulx_opt[0],uly_opt[0],uli,ulj); + // imgReader.geo2image(lrx_opt[0],lry_opt[0],lri,lrj); + // imgReader.getMinMax(static_cast<int>(uli),static_cast<int>(lri),static_cast<int>(ulj),static_cast<int>(lrj),band_opt[iband],minValue,maxValue); + // } + // else + // imgReader.getMinMax(minValue,maxValue,band_opt[iband],true); + // if(minmax_opt[0]) + // std::cout << "-min " << minValue << " -max " << maxValue << " "; + // else{ + // if(min_opt[0]) + // std::cout << "-min " << minValue << " "; + // if(max_opt[0]) + // std::cout << "-max " << maxValue << " "; + // } + // } + // } + // if(relative_opt[0]) + // hist_opt[0]=true; + // if(hist_opt[0]){ + // assert(band_opt[0]<imgReader.nrOfBand()); + // unsigned int nbin=(nbin_opt.size())? nbin_opt[0]:0; + // std::vector<unsigned long int> output; + // minValue=0; + // maxValue=0; + // //todo: optimize such that getMinMax is only called once... + // imgReader.getMinMax(minValue,maxValue,band_opt[0]); + + // if(src_min_opt.size()) + // minValue=src_min_opt[0]; + // if(src_max_opt.size()) + // maxValue=src_max_opt[0]; + // unsigned long int nsample=imgReader.getHistogram(output,minValue,maxValue,nbin,band_opt[0]); + // std::cout.precision(10); + // for(int bin=0;bin<nbin;++bin){ + // double binValue=0; + // if(nbin==maxValue-minValue+1) + // binValue=minValue+bin; + // else + // binValue=minValue+static_cast<double>(maxValue-minValue)*(bin+0.5)/nbin; + // std::cout << binValue << " "; + // if(relative_opt[0]) + // std::cout << 100.0*static_cast<double>(output[bin])/static_cast<double>(nsample) << std::endl; + // else + // std::cout << static_cast<double>(output[bin]) << std::endl; + // } + // } diff --git a/src/apps/pkstatascii.cc b/src/apps/pkstatascii.cc index 03c8fb9..5e5e0aa 100644 --- a/src/apps/pkstatascii.cc +++ b/src/apps/pkstatascii.cc @@ -123,6 +123,14 @@ int main(int argc, char *argv[]) exit(0);//help was invoked, stop processing } + if(src_min_opt.size()){ + while(src_min_opt.size()<col_opt.size()) + src_min_opt.push_back(src_min_opt[0]); + } + if(src_max_opt.size()){ + while(src_max_opt.size()<col_opt.size()) + src_max_opt.push_back(src_max_opt[0]); + } statfactory::StatFactory stat; if(rand_opt[0]>0){ gsl_rng* r=stat.getRandomGenerator(time(NULL)); @@ -287,8 +295,6 @@ int main(int argc, char *argv[]) //end test stat.distribution(dataVector[icol],dataVector[icol].begin(),dataVector[icol].end(),statVector[icol],nbin,minValue,maxValue,sigma); - //test - cout << "debug1" << endl; if(verbose_opt[0]) std::cout << "min and max values: " << minValue << ", " << maxValue << std::endl; } @@ -348,7 +354,7 @@ int main(int argc, char *argv[]) // if(kde_opt[0]>0) // sigma=kde_opt[0]; // else - sigma=1.06*sqrt(sqrt(stat.var(dataVector[0]))*sqrt(stat.var(dataVector[0])))*pow(dataVector[0].size(),-0.2); + sigma=1.06*sqrt(sqrt(stat.var(dataVector[0]))*sqrt(stat.var(dataVector[1])))*pow(dataVector[0].size(),-0.2); } assert(nbin); if(verbose_opt[0]){ diff --git a/src/apps/pkstatogr.cc b/src/apps/pkstatogr.cc index a31925e..c646d07 100644 --- a/src/apps/pkstatogr.cc +++ b/src/apps/pkstatogr.cc @@ -178,7 +178,7 @@ int main(int argc, char *argv[]) else binValue=minValue+static_cast<double>(maxValue-minValue)*(ibin+0.5)/nbin; cout << binValue << " "; - if(relative_opt[0]) + if(relative_opt[0]||kde_opt[0]) cout << 100.0*static_cast<double>(binData[ibin])/theData.size() << endl; else cout << binData[ibin] << endl; diff --git a/src/apps/pksvm.cc b/src/apps/pksvm.cc index edac81a..3126b17 100644 --- a/src/apps/pksvm.cc +++ b/src/apps/pksvm.cc @@ -820,51 +820,49 @@ int main(int argc, char *argv[]) } oldRowMask=rowMask; } - } - else - continue;//no coverage in this mask - short theMask=0; - for(short ivalue=0;ivalue<msknodata_opt.size();++ivalue){ - if(msknodata_opt[ivalue]>=0){//values set in msknodata_opt are invalid - if(lineMask[colMask]==msknodata_opt[ivalue]){ - theMask=lineMask[colMask]; - masked=true; + short theMask=0; + for(short ivalue=0;ivalue<msknodata_opt.size();++ivalue){ + if(msknodata_opt[ivalue]>=0){//values set in msknodata_opt are invalid + if(lineMask[colMask]==msknodata_opt[ivalue]){ + theMask=lineMask[colMask]; + masked=true; break; + } } - } - else{//only values set in msknodata_opt are valid - if(lineMask[colMask]!=-msknodata_opt[ivalue]){ - theMask=lineMask[colMask]; - masked=true; - } - else{ - masked=false; - break; + else{//only values set in msknodata_opt are valid + if(lineMask[colMask]!=-msknodata_opt[ivalue]){ + theMask=lineMask[colMask]; + masked=true; + } + else{ + masked=false; + break; + } } } + if(masked){ + if(classBag_opt.size()) + for(int ibag=0;ibag<nbag;++ibag) + classBag[ibag][icol]=theMask; + classOut[icol]=theMask; + continue; + } } - if(masked){ + bool valid=false; + for(int iband=0;iband<hpixel[icol].size();++iband){ + if(hpixel[icol][iband]){ + valid=true; + break; + } + } + if(!valid){ if(classBag_opt.size()) for(int ibag=0;ibag<nbag;++ibag) - classBag[ibag][icol]=theMask; - classOut[icol]=theMask; - continue; + classBag[ibag][icol]=nodata_opt[0]; + classOut[icol]=nodata_opt[0]; + continue;//next column } } - bool valid=false; - for(int iband=0;iband<hpixel[icol].size();++iband){ - if(hpixel[icol][iband]){ - valid=true; - break; - } - } - if(!valid){ - if(classBag_opt.size()) - for(int ibag=0;ibag<nbag;++ibag) - classBag[ibag][icol]=nodata_opt[0]; - classOut[icol]=nodata_opt[0]; - continue;//next column - } for(short iclass=0;iclass<nclass;++iclass) probOut[iclass][icol]=0; if(verbose_opt[0]>1) diff --git a/src/imageclasses/ImgReaderGdal.cc b/src/imageclasses/ImgReaderGdal.cc index 63d25bd..40ae36a 100644 --- a/src/imageclasses/ImgReaderGdal.cc +++ b/src/imageclasses/ImgReaderGdal.cc @@ -21,6 +21,7 @@ along with pktools. If not, see <http://www.gnu.org/licenses/>. #include <assert.h> #include <sstream> #include <iostream> +#include <gsl/gsl_cdf.h> ImgReaderGdal::ImgReaderGdal(void) : m_gds(NULL), m_ncol(0), m_nrow(0), m_nband(0) @@ -506,16 +507,34 @@ void ImgReaderGdal::getMinMax(double& minValue, double& maxValue, int band, bool } } -unsigned long int ImgReaderGdal::getHistogram(std::vector<unsigned long int>& histvector, double& min, double& max, unsigned int& nbin, int theBand) const{ +double ImgReaderGdal::getHistogram(std::vector<double>& histvector, double& min, double& max, unsigned int& nbin, int theBand, bool kde){ double minValue=0; double maxValue=0; - getMinMax(minValue,maxValue,theBand); - if(min<max&&min>minValue) + double meanValue=0; + double stdDev=0; + GDALProgressFunc pfnProgress; + void* pProgressData; + GDALRasterBand* rasterBand; + rasterBand=getRasterBand(theBand); + rasterBand->ComputeStatistics(0,&minValue,&maxValue,&meanValue,&stdDev,pfnProgress,pProgressData); + + if(min>=max) + getMinMax(minValue,maxValue,theBand); + else{ minValue=min; - if(min<max&&max<maxValue) maxValue=max; + } + // if(min<max&&min>minValue) + // minValue=min; + // if(min<max&&max<maxValue) + // maxValue=max; min=minValue; max=maxValue; + + double sigma=0; + if(kde) + sigma=1.06*stdDev*pow(getNvalid(theBand),-0.2); + double scale=0; if(maxValue>minValue){ if(nbin==0) @@ -526,6 +545,7 @@ unsigned long int ImgReaderGdal::getHistogram(std::vector<unsigned long int>& hi nbin=1; assert(nbin>0); histvector.resize(nbin); + double nvalid=0; unsigned long int nsample=0; unsigned long int ninvalid=0; std::vector<double> lineBuffer(nrOfCol()); @@ -542,10 +562,25 @@ unsigned long int ImgReaderGdal::getHistogram(std::vector<unsigned long int>& hi else if(nbin==1) ++histvector[0]; else{//scale to [0:nbin] - int theBin=static_cast<unsigned long int>(scale*(lineBuffer[icol]-minValue)); - assert(theBin>=0); - assert(theBin<nbin); - ++histvector[theBin]; + if(sigma>0){ + //create kde for Gaussian basis function + //todo: speed up by calculating first and last bin with non-zero contriubtion... + //todo: calculate real surface below pdf by using gsl_cdf_gaussian_P(x-mean+binsize,sigma)-gsl_cdf_gaussian_P(x-mean,sigma) + //hiero + for(int ibin=0;ibin<nbin;++ibin){ + double icenter=minValue+static_cast<double>(maxValue-minValue)*(ibin+0.5)/nbin; + double thePdf=gsl_ran_gaussian_pdf(lineBuffer[icol]-icenter, sigma); + histvector[ibin]+=thePdf; + nvalid+=thePdf; + } + } + else{ + int theBin=static_cast<unsigned long int>(scale*(lineBuffer[icol]-minValue)); + assert(theBin>=0); + assert(theBin<nbin); + ++histvector[theBin]; + ++nvalid; + } // else if(lineBuffer[icol]==maxValue) // ++histvector[nbin-1]; // else @@ -553,7 +588,7 @@ unsigned long int ImgReaderGdal::getHistogram(std::vector<unsigned long int>& hi } } } - unsigned long int nvalid=nrOfCol()*nrOfRow()-ninvalid; + // unsigned long int nvalid=nrOfCol()*nrOfRow()-ninvalid; return nvalid; } @@ -571,6 +606,26 @@ void ImgReaderGdal::getRange(std::vector<short>& range, int band) const sort(range.begin(),range.end()); } +unsigned long int ImgReaderGdal::getNvalid(int band) const +{ + unsigned long int nvalid=0; + if(m_noDataValues.size()){ + std::vector<short> lineBuffer(nrOfCol()); + for(int irow=0;irow<nrOfRow();++irow){ + readData(lineBuffer,GDT_Float64,irow,band); + for(int icol=0;icol<nrOfCol();++icol){ + if(isNoData(lineBuffer[icol])) + continue; + else + ++nvalid; + } + } + return nvalid; + } + else + return(nrOfCol()*nrOfRow()); +} + int ImgReaderGdal::getNoDataValues(std::vector<double>& noDataValues) const { if(m_noDataValues.size()){ diff --git a/src/imageclasses/ImgReaderGdal.h b/src/imageclasses/ImgReaderGdal.h index d287631..2d037e7 100644 --- a/src/imageclasses/ImgReaderGdal.h +++ b/src/imageclasses/ImgReaderGdal.h @@ -99,10 +99,11 @@ public: void getMinMax(int startCol, int endCol, int startRow, int endRow, int band, double& minValue, double& maxValue) const; void getMinMax(double& minValue, double& maxValue, int band=0, bool exhaustiveSearch=false) const; double getMin(int& col, int& row, int band=0) const; - unsigned long int getHistogram(std::vector<unsigned long int>& histvector, double& min, double& max,unsigned int& nbin, int theBand=0) const; + double getHistogram(std::vector<double>& histvector, double& min, double& max,unsigned int& nbin, int theBand=0, bool kde=false); double getMax(int& col, int& row, int band=0) const; void getRefPix(double& refX, double &refY, int band=0) const; void getRange(std::vector<short>& range, int Band=0) const; + unsigned long int getNvalid(int band) const; GDALDataType getDataType(int band=0) const; GDALRasterBand* getRasterBand(int band=0); GDALColorTable* getColorTable(int band=0) const; diff --git a/src/imageclasses/Makefile.am b/src/imageclasses/Makefile.am index 3f03c8b..216c7d2 100644 --- a/src/imageclasses/Makefile.am +++ b/src/imageclasses/Makefile.am @@ -15,7 +15,7 @@ libimageClasses_ladir = $(includedir)/pktools/imageclasses ## Instruct libtool to include ABI version information in the generated shared ## library file (.so). The library ABI version is defined in configure.ac, so ## that all version information is kept in one place. -libimageClasses_la_LDFLAGS = -version-info $(PKTOOLS_SO_VERSION) $(AM_LDFLAGS) +libimageClasses_la_LDFLAGS = -version-info $(PKTOOLS_SO_VERSION) $(AM_LDFLAGS) -lgsl # the list of header files that belong to the library (to be installed later) libimageClasses_la_HEADERS = ImgReaderGdal.h ImgReaderOgr.h ImgWriterGdal.h ImgWriterOgr.h diff --git a/src/imageclasses/Makefile.in b/src/imageclasses/Makefile.in index 6e446c5..1a61779 100644 --- a/src/imageclasses/Makefile.in +++ b/src/imageclasses/Makefile.in @@ -368,7 +368,7 @@ lib_LTLIBRARIES = libimageClasses.la # where to install the headers on the system libimageClasses_ladir = $(includedir)/pktools/imageclasses -libimageClasses_la_LDFLAGS = -version-info $(PKTOOLS_SO_VERSION) $(AM_LDFLAGS) +libimageClasses_la_LDFLAGS = -version-info $(PKTOOLS_SO_VERSION) $(AM_LDFLAGS) -lgsl # the list of header files that belong to the library (to be installed later) libimageClasses_la_HEADERS = ImgReaderGdal.h ImgReaderOgr.h ImgWriterGdal.h ImgWriterOgr.h -- Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/pkg-grass/pktools.git _______________________________________________ Pkg-grass-devel mailing list [email protected] http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-grass-devel
