[pktools] 01/10: Imported Upstream version 2.6.3
Sebastiaan Couwenberg
sebastic at moszumanska.debian.org
Sun Mar 8 15:57:05 UTC 2015
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 <sebastic at xs4all.nl>
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 <kempenep at gmail.com>.
#
@@ -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='kempenep at gmail.com'
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], [kempenep at gmail.com])
+AC_INIT([pktools], [2.6.3], [kempenep at gmail.com])
#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 at am__append_1 = pkann pkfsann pkregann
@USE_LAS_TRUE at am__append_2 = pklas2img
@USE_NLOPT_TRUE at am__append_3 = pkoptsvm
@@ -114,7 +112,7 @@ CONFIG_CLEAN_VPATH_FILES =
@USE_LAS_TRUE at am__EXEEXT_2 = pklas2img$(EXEEXT)
@USE_NLOPT_TRUE at 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 at 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 at ./$(DEPDIR)/pkregann-pkregann.Po at am__quote@
@AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/pksetmask.Po at am__quote@
@AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/pksieve.Po at am__quote@
+ at AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/pkstat.Po at am__quote@
@AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/pkstatascii.Po at am__quote@
@AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/pkstatogr.Po at am__quote@
@AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/pksvm.Po at 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
More information about the Pkg-grass-devel
mailing list