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

-- 
Alioth's /usr/local/bin/git-commit-notice on 
/srv/git.debian.org/git/pkg-grass/pktools.git

_______________________________________________
Pkg-grass-devel mailing list
Pkg-grass-devel@lists.alioth.debian.org
http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-grass-devel

Reply via email to