[pktools] 01/05: Imported Upstream version 2.6.2
Bas Couwenberg
sebastic at xs4all.nl
Wed Jan 28 23:26:53 UTC 2015
This is an automated email from the git hooks/post-receive script.
sebastic-guest pushed a commit to branch master
in repository pktools.
commit ac2bcd735505dc5de5efa66abdb438a1be4cf4d0
Author: Bas Couwenberg <sebastic at xs4all.nl>
Date: Wed Jan 28 09:21:07 2015 +0100
Imported Upstream version 2.6.2
---
ChangeLog | 10 +
configure | 20 +-
configure.ac | 2 +-
src/algorithms/Filter.cc | 181 ++++++++++++++++-
src/algorithms/Filter.h | 39 ++--
src/algorithms/Filter2d.h | 3 +-
src/algorithms/StatFactory.h | 58 ++++--
src/apps/pkann.cc | 29 ++-
src/apps/pkascii2img.cc | 5 +-
src/apps/pkascii2ogr.cc | 5 +-
src/apps/pkcomposite.cc | 42 +++-
src/apps/pkcreatect.cc | 10 +-
src/apps/pkcrop.cc | 42 +++-
src/apps/pkdiff.cc | 274 +++++++++++++-------------
src/apps/pkdsm2shadow.cc | 12 +-
src/apps/pkdumpimg.cc | 11 +-
src/apps/pkdumpogr.cc | 5 +-
src/apps/pkextract.cc | 57 ++++--
src/apps/pkfillnodata.cc | 10 +-
src/apps/pkfilter.cc | 158 +++++++++++----
src/apps/pkfilterascii.cc | 22 ++-
src/apps/pkfilterdem.cc | 20 +-
src/apps/pkfsann.cc | 30 ++-
src/apps/pkfssvm.cc | 36 +++-
src/apps/pkgetmask.cc | 18 +-
src/apps/pkinfo.cc | 3 +
src/apps/pkkalman.cc | 450 ++++++++++++++++++++++++++++---------------
src/apps/pklas2img.cc | 32 +--
src/apps/pkoptsvm.cc | 61 ++++--
src/apps/pkpolygonize.cc | 9 +-
src/apps/pkregann.cc | 19 +-
src/apps/pksetmask.cc | 16 +-
src/apps/pksieve.cc | 11 +-
src/apps/pkstatascii.cc | 32 ++-
src/apps/pkstatogr.cc | 7 +-
src/apps/pksvm.cc | 42 +++-
src/base/Optionpk.h | 26 ++-
37 files changed, 1311 insertions(+), 496 deletions(-)
diff --git a/ChangeLog b/ChangeLog
index bdc27ac..59a04c0 100755
--- a/ChangeLog
+++ b/ChangeLog
@@ -328,6 +328,7 @@ version 2.6.1
- pkfilter
Declare nodata option as double (see ticket #43500)
Support nodata values for filtering in spectral/temporal domain (see ticket #43713)
+ make use of memory buffer vsimem for morphological closing and opening instead of tmp file
- pkextract
extracting with absolute threshold (negative value) was selecting 1 extra sample unit
debug: removed (redundant?) SetGeometry after SetFrom when readFeature was of appropriate geometry type. This caused wrong polygons in output vectro dataset
@@ -336,6 +337,15 @@ version 2.6.1
- pkann
support for mask in geo coordinates (does not need to be in same dimensions as input raster dataset, only in same projection)
+version 2.6.2
+ - documentation
+ minor bug fixes, thanks to patches introduced by Bas Couwenberg
+ - pkdiff
+ support double data types for input and mask
+ support new option -reg to calculate regression between input and reference
+ 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)
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 86da08a..9106467 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.1.
+# Generated by GNU Autoconf 2.69 for pktools 2.6.2.
#
# Report bugs to <kempenep at gmail.com>.
#
@@ -590,8 +590,8 @@ MAKEFLAGS=
# Identity of this package.
PACKAGE_NAME='pktools'
PACKAGE_TARNAME='pktools'
-PACKAGE_VERSION='2.6.1'
-PACKAGE_STRING='pktools 2.6.1'
+PACKAGE_VERSION='2.6.2'
+PACKAGE_STRING='pktools 2.6.2'
PACKAGE_BUGREPORT='kempenep at gmail.com'
PACKAGE_URL=''
@@ -1366,7 +1366,7 @@ if test "$ac_init_help" = "long"; then
# Omit some internal or obsolete options to make the list less imposing.
# This message is too long to be a string in the A/UX 3.1 sh.
cat <<_ACEOF
-\`configure' configures pktools 2.6.1 to adapt to many kinds of systems.
+\`configure' configures pktools 2.6.2 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.1:";;
+ short | recursive ) echo "Configuration of pktools 2.6.2:";;
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.1
+pktools configure 2.6.2
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.1, which was
+It was created by pktools $as_me 2.6.2, 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.1'
+ VERSION='2.6.2'
cat >>confdefs.h <<_ACEOF
@@ -20172,7 +20172,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.1, which was
+This file was extended by pktools $as_me 2.6.2, which was
generated by GNU Autoconf 2.69. Invocation command line was
CONFIG_FILES = $CONFIG_FILES
@@ -20238,7 +20238,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.1
+pktools config.status 2.6.2
configured by $0, generated by GNU Autoconf 2.69,
with options \\"\$ac_cs_config\\"
diff --git a/configure.ac b/configure.ac
index 469ba42..631ff56 100644
--- a/configure.ac
+++ b/configure.ac
@@ -1,4 +1,4 @@
-AC_INIT([pktools], [2.6.1], [kempenep at gmail.com])
+AC_INIT([pktools], [2.6.2], [kempenep at gmail.com])
#AM_INIT_AUTOMAKE([-Wall -Werror foreign])
AM_INIT_AUTOMAKE([-Wall -Wno-extra-portability foreign])
#AM_INIT_AUTOMAKE([subdir-objects]) #not working due to bug in autoconf, see Debian list: Bug #752993)
diff --git a/src/algorithms/Filter.cc b/src/algorithms/Filter.cc
index f2b6aa0..150b59a 100644
--- a/src/algorithms/Filter.cc
+++ b/src/algorithms/Filter.cc
@@ -288,6 +288,39 @@ void filter::Filter::morphology(const ImgReaderGdal& input, ImgWriterGdal& outpu
}
}
+void filter::Filter::smoothNoData(const ImgReaderGdal& input, const std::string& interpolationType, ImgWriterGdal& output)
+{
+ Vector2d<double> lineInput(input.nrOfBand(),input.nrOfCol());
+ Vector2d<double> lineOutput(input.nrOfBand(),input.nrOfCol());
+ const char* pszMessage;
+ void* pProgressArg=NULL;
+ GDALProgressFunc pfnProgress=GDALTermProgress;
+ double progress=0;
+ pfnProgress(progress,pszMessage,pProgressArg);
+ for(int y=0;y<input.nrOfRow();++y){
+ for(int iband=0;iband<input.nrOfBand();++iband)
+ input.readData(lineInput[iband],GDT_Float64,y,iband);
+ vector<double> pixelInput(input.nrOfBand());
+ vector<double> pixelOutput(input.nrOfBand());
+ for(int x=0;x<input.nrOfCol();++x){
+ pixelInput=lineInput.selectCol(x);
+ smoothNoData(pixelInput,interpolationType,pixelOutput);
+ for(int iband=0;iband<input.nrOfBand();++iband)
+ lineOutput[iband][x]=pixelOutput[iband];
+ }
+ for(int iband=0;iband<input.nrOfBand();++iband){
+ try{
+ output.writeData(lineOutput[iband],GDT_Float64,y,iband);
+ }
+ catch(string errorstring){
+ cerr << errorstring << "in band " << iband << ", line " << y << endl;
+ }
+ }
+ progress=(1.0+y)/output.nrOfRow();
+ pfnProgress(progress,pszMessage,pProgressArg);
+ }
+}
+
void filter::Filter::smooth(const ImgReaderGdal& input, ImgWriterGdal& output, short dim)
{
assert(dim>0);
@@ -297,15 +330,6 @@ void filter::Filter::smooth(const ImgReaderGdal& input, ImgWriterGdal& output, s
filter(input,output);
}
-// void filter::Filter::smoothnodata(const ImgReaderGdal& input, ImgWriterGdal& output, short dim, short down, int offset)
-// {
-// assert(dim>0);
-// m_taps.resize(dim);
-// for(int itap=0;itap<dim;++itap)
-// m_taps[itap]=1.0/dim;
-// filter(input,output,down,offset);
-// }
-
void filter::Filter::filter(const ImgReaderGdal& input, ImgWriterGdal& output)
{
Vector2d<double> lineInput(input.nrOfBand(),input.nrOfCol());
@@ -432,6 +456,145 @@ void filter::Filter::filter(const ImgReaderGdal& input, ImgWriterGdal& output, c
}
}
+void filter::Filter::getSavGolayCoefficients(vector<double> &tapz, int np, int nl, int nr, int ld, int m) {
+ int j, k, imj, ipj, kk, mm;
+ double d, fac, sum;
+
+ // c.resize(nl+1+nr);
+ vector<double> tmpc(np);
+ if(np < nl + nr + 1 || nl < 0 || nr < 0 || ld > m || nl + nr < m) {
+ cerr << "bad args in savgol" << endl;
+ return;
+ }
+ vector<int> indx(m + 1, 0);
+ vector<double> a((m + 1) * (m + 1), 0.0);
+ vector<double> b(m + 1, 0.0);
+
+ for(ipj = 0; ipj <= (m << 1); ++ipj) {
+ sum = (ipj ? 0.0 : 1.0);
+ for(k = 1; k <= nr; ++k)
+ sum += (int)pow((double)k, (double)ipj);
+ for(k = 1; k <= nl; ++k)
+ sum += (int)pow((double) - k, (double)ipj);
+ mm = (ipj < 2 * m - ipj ? ipj : 2 * m - ipj);
+ for(imj = -mm; imj <= mm; imj += 2)
+ a[(ipj + imj) / 2 * (m + 1) + (ipj - imj) / 2] = sum;
+ }
+ ludcmp(a, indx, d);
+
+ for(j = 0; j < m + 1; ++j)
+ b[j] = 0.0;
+ b[ld] = 1.0;
+
+ lubksb(a, indx, b);
+ // for(kk = 0; kk < np; ++kk)
+ // c[kk] = 0.0;
+ for(k = -nl; k <= nr; ++k) {
+ // for(k = -nl; k < nr; ++k) {
+ sum = b[0];
+ fac = 1.0;
+ for(mm = 1; mm <= m; ++mm)
+ sum += b[mm] * (fac *= k);
+ // store in wrap=around order
+ kk = (np - k) % np;
+ //re-order c as I need for taps
+ // kk=k+nl;
+ tmpc[kk] = sum;
+ }
+ tapz.resize(nl+1+nr);
+ // for(k=0;k<nl+1+nr)
+ tapz[tapz.size()/2]=tmpc[0];
+ //past data points
+ for(k=1;k<=tapz.size()/2;++k)
+ tapz[tapz.size()/2-k]=tmpc[k];
+ //future data points
+ for(k=1;k<=tapz.size()/2;++k)
+ tapz[tapz.size()/2+k]=tmpc[np-k];
+}
+
+void filter::Filter::ludcmp(vector<double> &a, vector<int> &indx, double &d) {
+ const double TINY = 1.0e-20;
+ int i, imax = -1, j, k;
+ double big, dum, sum, temp;
+
+ int n = indx.size();
+ vector<double> vv(n, 0.0);
+
+ d = 1.0;
+ for(i = 0; i < n; ++i) {
+ big = 0.0;
+ for(j = 0; j < n; ++j)
+ if((temp = fabs(a[i * n + j])) > big) big = temp;
+
+ if(big == 0.0) {
+ cerr << "Singular matrix in routine ludcmp" << endl;
+ return;
+ }
+ vv[i] = 1. / big;
+ }
+
+ for(j = 0; j < n; ++j) {
+ for(i = 0; i < j; ++i) {
+ sum = a[i * n + j];
+ for(k = 0; k < i; ++k)
+ sum -= a[i * n + k] * a[k * n + j];
+ a[i * n + j] = sum;
+ }
+ big = 0.0;
+ for(i = j; i < n; ++i) {
+ sum = a[i * n + j];
+ for(k = 0; k < j; ++k)
+ sum -= a[i * n + k] * a[k * n + j];
+ a[i * n + j] = sum;
+ if((dum = vv[i] * fabs(sum)) >= big) {
+ big = dum;
+ imax = i;
+ }
+ }
+
+ if(j != imax) {
+ for(k = 0; k < n; ++k) {
+ dum = a[imax * n + k];
+ a[imax * n + k] = a[j * n + k];
+ a[j * n + k] = dum;
+ }
+ d = -d;
+ vv[imax] = vv[j];
+ }
+ indx[j] = imax;
+ if(a[j * n + j] == 0.0) a[j * n + j] = TINY;
+ if(j != n - 1) {
+ dum = 1. / a[j * n + j];
+ for(i = j + 1; i < n; ++i)
+ a[i * n + j] *= dum;
+ }
+ }
+}
+
+void filter::Filter::lubksb(vector<double> &a, vector<int> &indx, vector<double> &b) {
+ int i, ii = 0, ip, j;
+ double sum;
+ int n = indx.size();
+
+ for(i = 0; i < n; ++i) {
+ ip = indx[i];
+ sum = b[ip];
+ b[ip] = b[i];
+ if(ii != 0)
+ for(j = ii - 1; j < i; ++j)
+ sum -= a[i * n + j] * b[j];
+ else if(sum != 0.0)
+ ii = i + 1;
+ b[i] = sum;
+ }
+ for(i = n - 1; i >= 0; --i) {
+ sum = b[i];
+ for(j = i + 1; j < n; ++j)
+ sum -= a[i * n + j] * b[j];
+ b[i] = sum / a[i * n + i];
+ }
+}
+
double filter::Filter::getCentreWavelength(const std::vector<double> &wavelengthIn, const Vector2d<double>& srf, const std::string& interpolationType, double delta, bool verbose)
{
assert(srf.size()==2);//[0]: wavelength, [1]: response function
diff --git a/src/algorithms/Filter.h b/src/algorithms/Filter.h
index a75626c..171ce9d 100644
--- a/src/algorithms/Filter.h
+++ b/src/algorithms/Filter.h
@@ -33,9 +33,9 @@ extern "C" {
namespace filter
{
- enum FILTER_TYPE { median=0, var=1 , min=2, max=3, sum=4, mean=5, minmax=6, dilate=7, erode=8, close=9, open=10, homog=11, sobelx=12, sobely=13, sobelxy=14, sobelyx=-14, smooth=15, density=16, mode=17, mixed=18, smoothnodata=19, threshold=20, ismin=21, ismax=22, heterog=23, order=24, stdev=25, dwt=26, dwti=27, dwt_cut=28, dwt_cut_from=29};
+ enum FILTER_TYPE { median=0, var=1 , min=2, max=3, sum=4, mean=5, minmax=6, dilate=7, erode=8, close=9, open=10, homog=11, sobelx=12, sobely=13, sobelxy=14, sobelyx=-14, smooth=15, density=16, mode=17, mixed=18, smoothnodata=19, threshold=20, ismin=21, ismax=22, heterog=23, order=24, stdev=25, dwt=26, dwti=27, dwt_cut=28, dwt_cut_from=29, savgolay=30};
- enum PADDING { symmetric=0, replicate=1, circular=2, constant=3};
+ enum PADDING { symmetric=0, replicate=1, circular=2, zero=3};
class Filter
{
@@ -69,6 +69,7 @@ public:
template<class T> void filter(const std::vector<T>& input, std::vector<T>& output);
template<class T> void filter(const std::vector<T>& input, std::vector<T>& output, const std::string& method, int dim);
template<class T> void smooth(const std::vector<T>& input, std::vector<T>& output, short dim);
+ template<class T> void smoothNoData(const std::vector<T>& input, const std::string& interpolationType, std::vector<T>& output);
template<class T> void filter(T* input, int inputSize, std::vector<T>& output);
template<class T> void smooth(T* input, int inputSize, std::vector<T>& output, short dim);
//template<class T> void morphology(const std::vector<T>& input, std::vector<T>& output, const std::string& method, int dim, bool verbose=false);
@@ -76,7 +77,12 @@ public:
void filter(const ImgReaderGdal& input, ImgWriterGdal& output);
void stat(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method);
void filter(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, int dim);
+ void getSavGolayCoefficients(std::vector<double> &c, int np, int nl, int nr, int ld, int m);
+ void ludcmp(std::vector<double> &a, std::vector<int> &indx, double &d);
+ void lubksb(std::vector<double> &a, std::vector<int> &indx, std::vector<double> &b);
+ /* void savgolay(const ImgReaderGdal& input, ImgWriterGdal& output, int np, int nl, int nr, int m); */
void smooth(const ImgReaderGdal& input, ImgWriterGdal& output, short dim);
+ void smoothNoData(const ImgReaderGdal& input, const std::string& interpolationType, ImgWriterGdal& output);
double getCentreWavelength(const std::vector<double> &wavelengthIn, const Vector2d<double>& srf, const std::string& interpolationType, double delta=1.0, bool verbose=false);
template<class T> double applySrf(const std::vector<double> &wavelengthIn, const std::vector<T>& input, const Vector2d<double>& srf, const std::string& interpolationType, T& output, double delta=1.0, bool normalize=false, bool verbose=false);
template<class T> double applySrf(const std::vector<double> &wavelengthIn, const Vector2d<T>& input, const Vector2d<double>& srf, const std::string& interpolationType, std::vector<T>& output, double delta=1.0, bool normalize=false, int down=1, bool transposeInput=false, bool verbose=false);
@@ -126,12 +132,13 @@ private:
m_filterMap["heterog"]=filter::heterog;
m_filterMap["order"]=filter::order;
m_filterMap["median"]=filter::median;
+ m_filterMap["savgolay"]=filter::savgolay;
}
static PADDING getPadding(const std::string& padString){
std::map<std::string, PADDING> padMap;
- padMap["constant"]=filter::constant;
+ padMap["zero"]=filter::zero;
padMap["symmetric"]=filter::symmetric;
padMap["replicate"]=filter::replicate;
padMap["circular"]=filter::circular;
@@ -419,9 +426,19 @@ template<class T> void Filter::applyFwhm(const std::vector<double> &wavelengthIn
filter(input,output);
}
+ template<class T> void Filter::smoothNoData(const std::vector<T>& input, const std::string& interpolationType, std::vector<T>& output)
+{
+ statfactory::StatFactory stat;
+ stat.setNoDataValues(m_noDataValues);
+ std::vector<double> abscis(input.size());
+ for(int i=0;i<abscis.size();++i)
+ abscis[i]=i;
+ stat.interpolateNoData(abscis,input,interpolationType,output);
+ }
+
template<class T> void Filter::filter(const std::vector<T>& input, std::vector<T>& output)
{
- assert(input.size()>m_taps.size());
+ assert(input.size()>=m_taps.size());
output.resize(input.size());
int i=0;
//start: extend input by padding
@@ -440,7 +457,7 @@ template<class T> void Filter::filter(const std::vector<T>& input, std::vector<T
case(circular):
output[i]+=m_taps[m_taps.size()/2-t]*input[input.size()+i-t];
break;
- case(constant):
+ case(zero):
output[i]+=m_taps[m_taps.size()/2-t]*0;
break;
case(symmetric):
@@ -477,7 +494,7 @@ template<class T> void Filter::filter(const std::vector<T>& input, std::vector<T
case(circular):
output[i]+=m_taps[m_taps.size()/2+t]*input[t-1];
break;
- case(constant):
+ case(zero):
output[i]+=m_taps[m_taps.size()/2+t]*0;
break;
case(symmetric):
@@ -540,7 +557,7 @@ template<class T> void Filter::filter(const std::vector<T>& input, std::vector<T
case(circular):
theValue=input[input.size()+i-t];
break;
- case(constant):
+ case(zero):
theValue=0;
break;
case(symmetric):
@@ -668,7 +685,7 @@ template<class T> void Filter::filter(const std::vector<T>& input, std::vector<T
case(circular):
theValue=input[t-1];
break;
- case(constant):
+ case(zero):
theValue=0;
break;
case(symmetric):
@@ -728,7 +745,7 @@ template<class T> void Filter::filter(const std::vector<T>& input, std::vector<T
template<class T> void Filter::filter(T* input, int inputSize, std::vector<T>& output)
{
- assert(inputSize>m_taps.size());
+ assert(inputSize>=m_taps.size());
output.resize(inputSize);
int i=0;
@@ -749,7 +766,7 @@ template<class T> void Filter::filter(T* input, int inputSize, std::vector<T>& o
case(circular):
output[i]+=m_taps[m_taps.size()/2-t]*input[input.size()+i-t];
break;
- case(constant):
+ case(zero):
output[i]+=m_taps[m_taps.size()/2-t]*0;
break;
case(symmetric):
@@ -786,7 +803,7 @@ template<class T> void Filter::filter(T* input, int inputSize, std::vector<T>& o
case(circular):
output[i]+=m_taps[m_taps.size()/2+t]*input[t-1];
break;
- case(constant):
+ case(zero):
output[i]+=m_taps[m_taps.size()/2+t]*0;
break;
case(symmetric):
diff --git a/src/algorithms/Filter2d.h b/src/algorithms/Filter2d.h
index 221c1df..e03f0c1 100644
--- a/src/algorithms/Filter2d.h
+++ b/src/algorithms/Filter2d.h
@@ -58,7 +58,7 @@ extern "C" {
namespace filter2d
{
- enum FILTER_TYPE { median=100, var=101 , min=102, max=103, sum=104, mean=105, minmax=106, dilate=107, erode=108, close=109, open=110, homog=111, sobelx=112, sobely=113, sobelxy=114, sobelyx=115, smooth=116, density=117, mode=118, mixed=119, threshold=120, ismin=121, ismax=122, heterog=123, order=124, stdev=125, mrf=126, dwt=127, dwti=128, dwt_cut=129, scramble=130, shift=131, linearfeature=132, smoothnodata=133, countid=134, dwt_cut_from=135};
+ enum FILTER_TYPE { median=100, var=101 , min=102, max=103, sum=104, mean=105, minmax=106, dilate=107, erode=108, close=109, open=110, homog=111, sobelx=112, sobely=113, sobelxy=114, sobelyx=115, smooth=116, density=117, mode=118, mixed=119, threshold=120, ismin=121, ismax=122, heterog=123, order=124, stdev=125, mrf=126, dwt=127, dwti=128, dwt_cut=129, scramble=130, shift=131, linearfeature=132, smoothnodata=133, countid=134, dwt_cut_from=135, savgolay=136};
enum RESAMPLE { NEAR = 0, BILINEAR = 1, BICUBIC = 2 };//bicubic not supported yet...
@@ -167,6 +167,7 @@ private:
m_filterMap["shift"]=filter2d::shift;
m_filterMap["linearfeature"]=filter2d::linearfeature;
m_filterMap["countid"]=filter2d::countid;
+ m_filterMap["savgolay"]=filter2d::savgolay;
}
Vector2d<double> m_taps;
diff --git a/src/algorithms/StatFactory.h b/src/algorithms/StatFactory.h
index 4e575e3..b9511d8 100644
--- a/src/algorithms/StatFactory.h
+++ b/src/algorithms/StatFactory.h
@@ -88,8 +88,8 @@ public:
}
assert(spline);
};
- static void initSpline(gsl_spline *spline, const double *x, const double *y, int size){
- gsl_spline_init (spline, x, y, size);
+ static int initSpline(gsl_spline *spline, const double *x, const double *y, int size){
+ return gsl_spline_init (spline, x, y, size);
};
static double evalSpline(gsl_spline *spline, double x, gsl_interp_accel *acc){
return gsl_spline_eval (spline, x, acc);
@@ -179,6 +179,7 @@ public:
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;
+ template<class T> void interpolateNoData(const std::vector<double>& wavelengthIn, const std::vector<T>& input, const std::string& type, std::vector<T>& output, bool verbose=false) const;
template<class T> void interpolateUp(const std::vector<double>& wavelengthIn, const std::vector<T>& input, const std::vector<double>& wavelengthOut, const std::string& type, std::vector<T>& output, bool verbose=false) const;
template<class T> void interpolateUp(const std::vector<double>& wavelengthIn, const std::vector< std::vector<T> >& input, const std::vector<double>& wavelengthOut, const std::string& type, std::vector< std::vector<T> >& output, bool verbose=false) const;
// template<class T> void interpolateUp(const std::vector< std::vector<T> >& input, std::vector< std::vector<T> >& output, double start, double end, double step, const gsl_interp_type* type);
@@ -1050,7 +1051,42 @@ template<class T> double StatFactory::linear_regression_err(const std::vector<T>
//alternatively: use GNU scientific library:
// gsl_stats_correlation (const double data1[], const size_t stride1, const double data2[], const size_t stride2, const size_t n)
-//todo: nodata?
+template<class T> void StatFactory::interpolateNoData(const std::vector<double>& wavelengthIn, const std::vector<T>& input, const std::string& type, std::vector<T>& output, bool verbose) const{
+ assert(wavelengthIn.size());
+ std::vector<double> wavelengthOut=wavelengthIn;
+ std::vector<T> validIn=input;
+ assert(input.size()==wavelengthIn.size());
+ int nband=wavelengthIn.size();
+ output.clear();
+ //remove nodata from input and corresponding wavelengthIn
+ if(m_noDataValues.size()){
+ typename std::vector<T>::iterator itValue=validIn.begin();
+ typename std::vector<T>::iterator itWavelength=wavelengthOut.begin();
+ while(itValue!=validIn.end()&&itWavelength!=wavelengthOut.end()){
+ if(isNoData(*itValue)){
+ validIn.erase(itValue);
+ wavelengthOut.erase(itWavelength);
+ }
+ else{
+ ++itValue;
+ ++itWavelength;
+ }
+ }
+ if(validIn.size()>1){
+ try{
+ interpolateUp(wavelengthOut, validIn, wavelengthIn, type, output, verbose);
+ }
+ catch(...){
+ output=input;
+ }
+ }
+ else//we can not interpolate if no valid data
+ output=input;
+ }
+ else//no nodata values to interpolate
+ output=input;
+}
+
template<class T> void StatFactory::interpolateUp(const std::vector<double>& wavelengthIn, const std::vector<T>& input, const std::vector<double>& wavelengthOut, const std::string& type, std::vector<T>& output, bool verbose) const{
assert(wavelengthIn.size());
assert(input.size()==wavelengthIn.size());
@@ -1064,7 +1100,11 @@ template<class T> void StatFactory::interpolateUp(const std::vector<double>& wav
assert(spline);
assert(&(wavelengthIn[0]));
assert(&(input[0]));
- initSpline(spline,&(wavelengthIn[0]),&(input[0]),nband);
+ int status=initSpline(spline,&(wavelengthIn[0]),&(input[0]),nband);
+ if(status){
+ std::string errorString="Could not initialize spline";
+ throw(errorString);
+ }
for(int index=0;index<wavelengthOut.size();++index){
if(wavelengthOut[index]<*wavelengthIn.begin()){
output.push_back(*(input.begin()));
@@ -1074,16 +1114,6 @@ template<class T> void StatFactory::interpolateUp(const std::vector<double>& wav
output.push_back(input.back());
continue;
}
- // if(type=="linear"){
- // if(wavelengthOut[index]<wavelengthIn.back()){
- // output.push_back(*(input.begin()));
- // continue;
- // }
- // else if(wavelengthOut[index]>wavelengthIn.back()){
- // output.push_back(input.back());
- // continue;
- // }
- // }
double dout=evalSpline(spline,wavelengthOut[index],acc);
output.push_back(dout);
}
diff --git a/src/apps/pkann.cc b/src/apps/pkann.cc
index b8d2ec5..95cc3fb 100644
--- a/src/apps/pkann.cc
+++ b/src/apps/pkann.cc
@@ -56,8 +56,8 @@ int main(int argc, char *argv[])
Optionpk<unsigned short> cv_opt("cv", "cv", "n-fold cross validation mode",0);
Optionpk<unsigned int> nneuron_opt("nn", "nneuron", "number of neurons in hidden layers in neural network (multiple hidden layers are set by defining multiple number of neurons: -n 15 -n 1, default is one hidden layer with 5 neurons)", 5);
Optionpk<float> connection_opt("\0", "connection", "connection reate (default: 1.0 for a fully connected network)", 1.0);
- Optionpk<float> weights_opt("w", "weights", "weights for neural network. Apply to fully connected network only, starting from first input neuron to last output neuron, including the bias neurons (last neuron in each but last layer)", 0.0);
Optionpk<float> learning_opt("l", "learning", "learning rate (default: 0.7)", 0.7);
+ Optionpk<float> weights_opt("w", "weights", "weights for neural network. Apply to fully connected network only, starting from first input neuron to last output neuron, including the bias neurons (last neuron in each but last layer)", 0.0);
Optionpk<unsigned int> maxit_opt("\0", "maxit", "number of maximum iterations (epoch) (default: 500)", 500);
Optionpk<unsigned short> comb_opt("comb", "comb", "how to combine bootstrap aggregation classifiers (0: sum rule, 1: product rule, 2: max rule). Also used to aggregate classes with rc option. Default is sum rule (0)",0);
Optionpk<unsigned short> bag_opt("bag", "bag", "Number of bootstrap aggregations (default is no bagging: 1)", 1);
@@ -78,7 +78,27 @@ int main(int argc, char *argv[])
Optionpk<unsigned int> nactive_opt("na", "nactive", "number of active training points",1);
Optionpk<string> classname_opt("c", "class", "list of class names.");
Optionpk<short> classvalue_opt("r", "reclass", "list of class values (use same order as in class opt).");
- Optionpk<short> verbose_opt("v", "verbose", "set to: 0 (results only), 1 (confusion matrix), 2 (debug)",0);
+ Optionpk<short> verbose_opt("v", "verbose", "set to: 0 (results only), 1 (confusion matrix), 2 (debug)",0,2);
+
+ band_opt.setHide(1);
+ bstart_opt.setHide(1);
+ bend_opt.setHide(1);
+ balance_opt.setHide(1);
+ minSize_opt.setHide(1);
+ bag_opt.setHide(1);
+ bagSize_opt.setHide(1);
+ comb_opt.setHide(1);
+ classBag_opt.setHide(1);
+ minSize_opt.setHide(1);
+ prob_opt.setHide(1);
+ priorimg_opt.setHide(1);
+ minSize_opt.setHide(1);
+ offset_opt.setHide(1);
+ scale_opt.setHide(1);
+ connection_opt.setHide(1);
+ weights_opt.setHide(1);
+ maxit_opt.setHide(1);
+ learning_opt.setHide(1);
bool doProcess;//stop process when program was invoked with help option (-h --help)
try{
@@ -129,7 +149,10 @@ int main(int argc, char *argv[])
exit(0);
}
if(!doProcess){
- std::cout << "short option -h shows basic options only, use long option --help to show all options" << std::endl;
+ cout << endl;
+ cout << "Usage: pkann -t training [-i input -o output] [-cv value]" << endl;
+ cout << endl;
+ cout << "short option -h shows basic options only, use long option --help to show all options" << endl;
exit(0);//help was invoked, stop processing
}
diff --git a/src/apps/pkascii2img.cc b/src/apps/pkascii2img.cc
index 26a6122..fbf24a6 100644
--- a/src/apps/pkascii2img.cc
+++ b/src/apps/pkascii2img.cc
@@ -39,7 +39,7 @@ int main(int argc, char *argv[])
Optionpk<string> projection_opt("a_srs", "a_srs", "Override the projection for the output file");
Optionpk<string> colorTable_opt("ct", "ct", "color table (file with 5 columns: id R G B ALFA (0: transparent, 255: solid)");
Optionpk<string> description_opt("d", "description", "Set image description");
- Optionpk<bool> verbose_opt("v", "verbose", "verbose", false);
+ Optionpk<bool> verbose_opt("v", "verbose", "verbose", false,2);
bool doProcess;//stop process when program was invoked with help option (-h --help)
try{
@@ -62,6 +62,9 @@ int main(int argc, char *argv[])
exit(0);
}
if(!doProcess){
+ cout << endl;
+ cout << "Usage: pkascii2img -i input.txt -o output" << 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
}
diff --git a/src/apps/pkascii2ogr.cc b/src/apps/pkascii2ogr.cc
index ddc648b..1cc32b3 100644
--- a/src/apps/pkascii2ogr.cc
+++ b/src/apps/pkascii2ogr.cc
@@ -37,7 +37,7 @@ int main(int argc, char *argv[])
Optionpk<string> ftype_opt("ot", "ot", "Field type (Real, Integer, String) for each of the fields as defined by name","Real");
Optionpk<string> projection_opt("a_srs", "a_srs", "Override the projection for the output file, use epsg:<code> or Wkt string", "epsg:4326");
Optionpk<char> fs_opt("fs","fs","field separator.",' ');
- Optionpk<int> verbose_opt("v", "verbose", "verbose (0)", 0);
+ Optionpk<int> verbose_opt("v", "verbose", "verbose (0)", 0,2);
bool doProcess;//stop process when program was invoked with help option (-h --help)
try{
@@ -58,6 +58,9 @@ int main(int argc, char *argv[])
exit(0);
}
if(!doProcess){
+ cout << endl;
+ cout << "Usage: pkascii2ogr -i input.txt -o output" << 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
}
diff --git a/src/apps/pkcomposite.cc b/src/apps/pkcomposite.cc
index 04e9897..2c4aeb0 100644
--- a/src/apps/pkcomposite.cc
+++ b/src/apps/pkcomposite.cc
@@ -29,7 +29,7 @@ along with pktools. If not, see <http://www.gnu.org/licenses/>.
#include "algorithms/StatFactory.h"
namespace crule{
- enum CRULE_TYPE {overwrite=0, maxndvi=1, maxband=2, minband=3, validband=4, mean=5, mode=6, median=7,sum=8};
+ enum CRULE_TYPE {overwrite=0, maxndvi=1, maxband=2, minband=3, validband=4, mean=5, mode=6, median=7,sum=8,minallbands=9,maxallbands=10};
}
using namespace std;
@@ -46,7 +46,7 @@ int main(int argc, char *argv[])
Optionpk<double> uly_opt("uly", "uly", "Upper left y value bounding box", 0.0);
Optionpk<double> lrx_opt("lrx", "lrx", "Lower right x value bounding box", 0.0);
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", "overwrite");
+ 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<int> bndnodata_opt("bndnodata", "bndnodata", "Bands in input image to check if pixel is valid (used for srcnodata, min and max options)", 0);
@@ -63,7 +63,14 @@ int main(int argc, char *argv[])
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)");
Optionpk<string> description_opt("d", "description", "Set image description");
- Optionpk<bool> verbose_opt("v", "verbose", "verbose", false);
+ Optionpk<bool> verbose_opt("v", "verbose", "verbose", false,2);
+
+ file_opt.setHide(1);
+ weight_opt.setHide(1);
+ class_opt.setHide(1);
+ colorTable_opt.setHide(1);
+ description_opt.setHide(1);
+ verbose_opt.setHide(1);
bool doProcess;//stop process when program was invoked with help option (-h --help)
try{
@@ -101,6 +108,9 @@ int main(int argc, char *argv[])
exit(0);
}
if(!doProcess){
+ cout << endl;
+ cout << "Usage: pkcomposite -i input [-i input]* -o output" << 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
}
@@ -118,6 +128,8 @@ int main(int argc, char *argv[])
cruleMap["mode"]=crule::mode;
cruleMap["median"]=crule::median;
cruleMap["sum"]=crule::sum;
+ cruleMap["maxallbands"]=crule::maxallbands;
+ cruleMap["minallbands"]=crule::minallbands;
while(srcnodata_opt.size()<bndnodata_opt.size())
srcnodata_opt.push_back(srcnodata_opt[0]);
@@ -264,6 +276,12 @@ int main(int argc, char *argv[])
case(crule::sum):
cout << "Composite rule: sum" << endl;
break;
+ case(crule::minallbands):
+ cout << "Composite rule: minallbands" << endl;
+ break;
+ case(crule::maxallbands):
+ cout << "Composite rule: maxallbands" << endl;
+ break;
}
}
if(band_opt.size()){
@@ -459,7 +477,7 @@ int main(int argc, char *argv[])
Vector2d< vector<double> > storeBuffer;
vector<bool> writeValid(ncol);
- if(cruleMap[crule_opt[0]]==crule::mean||cruleMap[crule_opt[0]]==crule::median||cruleMap[crule_opt[0]]==crule::sum)//mean, median or (weighted) sum value
+ if(cruleMap[crule_opt[0]]==crule::mean||cruleMap[crule_opt[0]]==crule::median||cruleMap[crule_opt[0]]==crule::sum||cruleMap[crule_opt[0]]==crule::minallbands||cruleMap[crule_opt[0]]==crule::maxallbands)//mean, median, (weighted) sum value, minallbands, maxallbands
storeBuffer.resize(nband,ncol);
for(int icol=0;icol<imgWriter.nrOfCol();++icol){
writeValid[icol]=false;
@@ -720,6 +738,8 @@ int main(int argc, char *argv[])
case(crule::mean)://mean value
case(crule::median)://median value
case(crule::sum)://sum value
+ case(crule::minallbands)://minimum for each and every band
+ case(crule::maxallbands)://maximum for each and every band
switch(theResample){
case(BILINEAR):
lowerCol=readCol-0.5;
@@ -781,7 +801,7 @@ int main(int argc, char *argv[])
++fileBuffer[ib];
break;
}
- }
+ }
else{
writeValid[ib]=true;//readValid was true
int iband=0;
@@ -789,6 +809,8 @@ int main(int argc, char *argv[])
case(crule::mean):
case(crule::median):
case(crule::sum):
+ case(crule::minallbands):
+ case(crule::maxallbands):
switch(theResample){
case(BILINEAR):
lowerCol=readCol-0.5;
@@ -932,6 +954,16 @@ int main(int argc, char *argv[])
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(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(storeBuffer[bands[iband]][icol].size())
+ writeBuffer[iband][icol]=stat.mymax(storeBuffer[bands[iband]][icol]);
+ break;
default:
break;
}
diff --git a/src/apps/pkcreatect.cc b/src/apps/pkcreatect.cc
index 295cb6f..0f8ae51 100644
--- a/src/apps/pkcreatect.cc
+++ b/src/apps/pkcreatect.cc
@@ -41,8 +41,11 @@ 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", "GTiff");
Optionpk<string> option_opt("co", "co", "Creation option for output file. Multiple options can be specified.");
Optionpk<string> description_opt("d", "description", "Set image description");
- Optionpk<bool> verbose_opt("v", "verbose", "verbose", false);
-
+ Optionpk<bool> verbose_opt("v", "verbose", "verbose", false,2);
+
+ legend_opt.setHide(1);
+ dim_opt.setHide(1);
+
bool doProcess;//stop process when program was invoked with help option (-h --help)
try{
doProcess=input_opt.retrieveOption(argc,argv);
@@ -63,6 +66,9 @@ int main(int argc,char **argv) {
exit(0);
}
if(!doProcess){
+ cout << endl;
+ cout << "Usage: pkcreatect -i input.txt -o output [-ct colortable | -min value -max value]" << 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
}
diff --git a/src/apps/pkcrop.cc b/src/apps/pkcrop.cc
index d926b91..6d6c1e9 100644
--- a/src/apps/pkcrop.cc
+++ b/src/apps/pkcrop.cc
@@ -61,37 +61,51 @@ int main(int argc, char *argv[])
Optionpk<double> nodata_opt("nodata", "nodata", "Nodata value to put in image if out of bounds.");
Optionpk<string> resample_opt("r", "resampling-method", "Resampling method (near: nearest neighbor, bilinear: bi-linear interpolation).", "near");
Optionpk<string> description_opt("d", "description", "Set image description");
- Optionpk<short> verbose_opt("v", "verbose", "verbose", 0);
+ Optionpk<short> verbose_opt("v", "verbose", "verbose", 0,2);
+ extent_opt.setHide(1);
+ mask_opt.setHide(1);
+ option_opt.setHide(1);
+ cx_opt.setHide(1);
+ cy_opt.setHide(1);
+ nx_opt.setHide(1);
+ ny_opt.setHide(1);
+ ns_opt.setHide(1);
+ nl_opt.setHide(1);
+ scale_opt.setHide(1);
+ offset_opt.setHide(1);
+ nodata_opt.setHide(1);
+ description_opt.setHide(1);
+
bool doProcess;//stop process when program was invoked with help option (-h --help)
try{
doProcess=input_opt.retrieveOption(argc,argv);
output_opt.retrieveOption(argc,argv);
projection_opt.retrieveOption(argc,argv);
- extent_opt.retrieveOption(argc,argv);
- mask_opt.retrieveOption(argc,argv);
ulx_opt.retrieveOption(argc,argv);
uly_opt.retrieveOption(argc,argv);
lrx_opt.retrieveOption(argc,argv);
lry_opt.retrieveOption(argc,argv);
band_opt.retrieveOption(argc,argv);
autoscale_opt.retrieveOption(argc,argv);
- scale_opt.retrieveOption(argc,argv);
- offset_opt.retrieveOption(argc,argv);
otype_opt.retrieveOption(argc,argv);
oformat_opt.retrieveOption(argc,argv);
- option_opt.retrieveOption(argc,argv);
colorTable_opt.retrieveOption(argc,argv);
dx_opt.retrieveOption(argc,argv);
dy_opt.retrieveOption(argc,argv);
+ resample_opt.retrieveOption(argc,argv);
+ extent_opt.retrieveOption(argc,argv);
+ mask_opt.retrieveOption(argc,argv);
+ option_opt.retrieveOption(argc,argv);
cx_opt.retrieveOption(argc,argv);
cy_opt.retrieveOption(argc,argv);
nx_opt.retrieveOption(argc,argv);
ny_opt.retrieveOption(argc,argv);
ns_opt.retrieveOption(argc,argv);
nl_opt.retrieveOption(argc,argv);
+ scale_opt.retrieveOption(argc,argv);
+ offset_opt.retrieveOption(argc,argv);
nodata_opt.retrieveOption(argc,argv);
- resample_opt.retrieveOption(argc,argv);
description_opt.retrieveOption(argc,argv);
verbose_opt.retrieveOption(argc,argv);
}
@@ -99,7 +113,14 @@ int main(int argc, char *argv[])
std::cout << predefinedString << std::endl;
exit(0);
}
+ //test
+ if(verbose_opt[0])
+ cout << setprecision(12) << "--ulx=" << ulx_opt[0] << " --uly=" << uly_opt[0] << " --lrx=" << lrx_opt[0] << " --lry=" << lry_opt[0] << endl;
+
if(!doProcess){
+ cout << endl;
+ cout << "Usage: pkcrop -i input -o output" << 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
}
@@ -190,6 +211,10 @@ int main(int argc, char *argv[])
double croplry=lry_opt[0];
//get bounding box from extentReader if defined
ImgReaderOgr extentReader;
+ //test
+ if(verbose_opt[0])
+ cout << "--ulx=" << ulx_opt[0] << " --uly=" << uly_opt[0] << " --lrx=" << lrx_opt[0] << " --lry=" << lry_opt[0] << endl;
+
if(extent_opt.size()){
for(int iextent=0;iextent<extent_opt.size();++iextent){
extentReader.open(extent_opt[iextent]);
@@ -214,6 +239,9 @@ int main(int argc, char *argv[])
lrx_opt[0]=cx_opt[0]+ns_opt[0]*dx/2.0;
lry_opt[0]=(isGeoRef) ? cy_opt[0]-nl_opt[0]*dy/2.0 : cy_opt[0]+nl_opt[0]*dy/2.0;
}
+ //test
+ if(verbose_opt[0])
+ cout << "--ulx=" << ulx_opt[0] << " --uly=" << uly_opt[0] << " --lrx=" << lrx_opt[0] << " --lry=" << lry_opt[0] << endl;
if(ulx_opt[0]<cropulx)
cropulx=ulx_opt[0];
if(uly_opt[0]>cropuly)
diff --git a/src/apps/pkdiff.cc b/src/apps/pkdiff.cc
index bb370be..375c591 100644
--- a/src/apps/pkdiff.cc
+++ b/src/apps/pkdiff.cc
@@ -32,26 +32,37 @@ int main(int argc, char *argv[])
Optionpk<string> input_opt("i", "input", "Input raster dataset.");
Optionpk<string> reference_opt("ref", "reference", "Reference (raster or vector) dataset");
Optionpk<string> layer_opt("ln", "ln", "Layer name(s) in sample. Leave empty to select all (for vector reference datasets only)");
- Optionpk<string> output_opt("o", "output", "Output dataset (optional)");
- Optionpk<string> ogrformat_opt("f", "f", "OGR format for output vector (for vector reference datasets only)","SQLite");
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<int> masknodata_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> valueE_opt("\0", "correct", "Value for correct pixels", 0,2);
- Optionpk<short> valueO_opt("\0", "omission", "Value for omission errors: input label > reference label", 1,2);
- Optionpk<short> valueC_opt("\0", "commission", "Value for commission errors: input label < reference label", 2,1);
+ 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<short> band_opt("b", "band", "Input raster band", 0);
+ 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);
Optionpk<bool> confusion_opt("cm", "confusion", "Create confusion matrix (to std out)", false);
Optionpk<string> labelref_opt("lr", "lref", "Attribute name of the reference label (for vector reference datasets only)", "label");
+ Optionpk<string> classname_opt("c", "class", "List of class names.");
+ Optionpk<short> classvalue_opt("r", "reclass", "List of class values (use same order as in classname option).");
+ Optionpk<string> output_opt("o", "output", "Output dataset (optional)");
+ Optionpk<string> ogrformat_opt("f", "f", "OGR format for output vector (for vector reference datasets only)","SQLite");
Optionpk<string> labelclass_opt("lc", "lclass", "Attribute name of the classified label (for vector reference datasets only)", "class");
Optionpk<short> boundary_opt("bnd", "boundary", "Boundary for selecting the sample (for vector reference datasets only)", 1,1);
Optionpk<bool> homogeneous_opt("hom", "homogeneous", "Only take regions with homogeneous boundary into account (for reference datasets only)", false,1);
Optionpk<bool> disc_opt("circ", "circular", "Use circular boundary (for vector reference datasets only)", false,1);
- Optionpk<string> classname_opt("c", "class", "List of class names.");
- Optionpk<short> classvalue_opt("r", "reclass", "List of class values (use same order as in classname option).");
Optionpk<string> colorTable_opt("ct", "ct", "Color table in ASCII format having 5 columns: id R G B ALFA (0: transparent, 255: solid).");
Optionpk<string> option_opt("co", "co", "Creation option for output file. Multiple options can be specified.");
- Optionpk<short> verbose_opt("v", "verbose", "Verbose level", 0);
+ Optionpk<short> valueE_opt("\0", "correct", "Value for correct pixels", 0,2);
+ Optionpk<short> valueO_opt("\0", "omission", "Value for omission errors: input label > reference label", 1,2);
+ Optionpk<short> valueC_opt("\0", "commission", "Value for commission errors: input label < reference label", 2,1);
+ Optionpk<short> verbose_opt("v", "verbose", "Verbose level", 0,2);
+
+ output_opt.setHide(1);
+ ogrformat_opt.setHide(1);
+ labelclass_opt.setHide(1);
+ boundary_opt.setHide(1);
+ homogeneous_opt.setHide(1);
+ disc_opt.setHide(1);
+ colorTable_opt.setHide(1);
+ option_opt.setHide(1);
bool doProcess;//stop process when program was invoked with help option (-h --help)
try{
@@ -59,25 +70,27 @@ int main(int argc, char *argv[])
reference_opt.retrieveOption(argc,argv);
layer_opt.retrieveOption(argc,argv);
band_opt.retrieveOption(argc,argv);
+ rmse_opt.retrieveOption(argc,argv);
+ regression_opt.retrieveOption(argc,argv);
confusion_opt.retrieveOption(argc,argv);
labelref_opt.retrieveOption(argc,argv);
classname_opt.retrieveOption(argc,argv);
classvalue_opt.retrieveOption(argc,argv);
nodata_opt.retrieveOption(argc,argv);
mask_opt.retrieveOption(argc,argv);
- masknodata_opt.retrieveOption(argc,argv);
+ msknodata_opt.retrieveOption(argc,argv);
output_opt.retrieveOption(argc,argv);
ogrformat_opt.retrieveOption(argc,argv);
labelclass_opt.retrieveOption(argc,argv);
- valueE_opt.retrieveOption(argc,argv);
- valueO_opt.retrieveOption(argc,argv);
- valueC_opt.retrieveOption(argc,argv);
boundary_opt.retrieveOption(argc,argv);
homogeneous_opt.retrieveOption(argc,argv);
disc_opt.retrieveOption(argc,argv);
colorTable_opt.retrieveOption(argc,argv);
option_opt.retrieveOption(argc,argv);
// class_opt.retrieveOption(argc,argv);
+ valueE_opt.retrieveOption(argc,argv);
+ valueO_opt.retrieveOption(argc,argv);
+ valueC_opt.retrieveOption(argc,argv);
verbose_opt.retrieveOption(argc,argv);
}
catch(string predefinedString){
@@ -85,6 +98,9 @@ int main(int argc, char *argv[])
exit(0);
}
if(!doProcess){
+ cout << endl;
+ cout << "Usage: pkdiff -i input -ref reference" << 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
}
@@ -108,6 +124,11 @@ int main(int argc, char *argv[])
exit(0);
}
+ //band_opt[0] is for input
+ //band_opt[1] is for reference
+ if(band_opt.size()<2)
+ band_opt.push_back(band_opt[0]);
+
if(mask_opt.size())
while(mask_opt.size()<input_opt.size())
mask_opt.push_back(mask_opt[0]);
@@ -118,6 +139,14 @@ int main(int argc, char *argv[])
map<string,short> classValueMap;
vector<std::string> nameVector(255);//the inverse of the classValueMap
vector<string> classNames;
+
+ unsigned int ntotalValidation=0;
+ unsigned int nflagged=0;
+ Vector2d<int> resultClass;
+ vector<float> user;
+ vector<float> producer;
+ vector<unsigned int> nvalidation;
+
if(confusion_opt[0]){
// if(class_opt.size()>1)
// inputRange=class_opt;
@@ -168,16 +197,6 @@ int main(int argc, char *argv[])
for(int iclass=0;iclass<cm.nClasses();++iclass)
cout << iclass << " " << cm.getClass(iclass) << endl;
}
- }
-
- unsigned int ntotalValidation=0;
- unsigned int nflagged=0;
- Vector2d<int> resultClass(nclass,nclass);
- vector<float> user(nclass);
- vector<float> producer(nclass);
- vector<unsigned int> nvalidation(nclass);
-
- if(confusion_opt[0]){
resultClass.resize(nclass,nclass);
user.resize(nclass);
producer.resize(nclass);
@@ -191,8 +210,8 @@ int main(int argc, char *argv[])
}
bool isDifferent=false;
-
bool refIsRaster=false;
+
ImgReaderOgr referenceReaderOgr;
try{
referenceReaderOgr.open(reference_opt[0]);
@@ -403,7 +422,7 @@ int main(int argc, char *argv[])
continue;
if(verbose_opt[0])
cout << setprecision(12) << "reading image value at x,y " << x << "," << y << " (" << i << "," << j << "), ";
- inputReader.readData(inputValue,GDT_Int16,i,j,band_opt[0]);
+ inputReader.readData(inputValue,GDT_Float64,i,j,band_opt[0]);
inputValues.push_back(inputValue);
if(inputValues.back()!=*(inputValues.begin()))
isHomogeneous=false;
@@ -416,18 +435,18 @@ int main(int argc, char *argv[])
break;
}
}
- maskFlagged=false;//(masknodata_opt[ivalue]>=0)?false:true;
+ maskFlagged=false;//(msknodata_opt[ivalue]>=0)?false:true;
if(mask_opt.size()){
- maskReader.readData(maskValue,GDT_Int16,i,j,band_opt[0]);
- for(int ivalue=0;ivalue<masknodata_opt.size();++ivalue){
- if(masknodata_opt[ivalue]>=0){//values set in masknodata_opt are invalid
- if(maskValue==masknodata_opt[ivalue]){
+ maskReader.readData(maskValue,GDT_Float64,i,j,0);
+ for(int ivalue=0;ivalue<msknodata_opt.size();++ivalue){
+ if(msknodata_opt[ivalue]>=0){//values set in msknodata_opt are invalid
+ if(maskValue==msknodata_opt[ivalue]){
maskFlagged=true;
break;
}
}
- else{//only values set in masknodata_opt are valid
- if(maskValue!=-masknodata_opt[ivalue])
+ else{//only values set in msknodata_opt are valid
+ if(maskValue!=-msknodata_opt[ivalue])
maskFlagged=true;
else{
maskFlagged=false;
@@ -571,11 +590,11 @@ int main(int argc, char *argv[])
}//next reference
}//next input
pfnProgress(1.0,pszMessage,pProgressArg);
- }//reference is OGR
- else{
+ }//reference is OGR vector
+ else{//reference is GDAL raster
ImgWriterGdal gdalWriter;
try{
- inputReader.open(input_opt[0]);//,imagicX_opt[0],imagicY_opt[0]);
+ inputReader.open(input_opt[0]);
if(mask_opt.size())
maskReader.open(mask_opt[0]);
if(output_opt.size()){
@@ -607,15 +626,18 @@ int main(int argc, char *argv[])
exit(2);
}
//todo: support different data types!
- vector<short> lineInput(inputReader.nrOfCol());
- vector<short> lineMask(maskReader.nrOfCol());
- vector<short> lineOutput;
+ vector<double> lineInput(inputReader.nrOfCol());
+ vector<double> lineMask(maskReader.nrOfCol());
+ vector<double> lineOutput;
+ vector<double> bufferInput;//for regression
+ vector<double> bufferReference;//for regression
if(output_opt.size())
lineOutput.resize(inputReader.nrOfCol());
int irow=0;
int icol=0;
double oldreferencerow=-1;
+ double oldmaskrow=-1;
ImgReaderGdal referenceReaderGdal;
try{
referenceReaderGdal.open(reference_opt[0]);//,rmagicX_opt[0],rmagicY_opt[0]);
@@ -627,11 +649,11 @@ int main(int argc, char *argv[])
if(inputReader.isGeoRef()){
assert(referenceReaderGdal.isGeoRef());
if(inputReader.getProjection()!=referenceReaderGdal.getProjection())
- cout << "projection of input image and reference image are different!" << endl;
+ cerr << "Warning: projection of input image and reference image are different" << endl;
}
- vector<short> lineReference(referenceReaderGdal.nrOfCol());
+ vector<double> lineReference(referenceReaderGdal.nrOfCol());
if(confusion_opt[0]){
- referenceReaderGdal.getRange(referenceRange,band_opt[0]);
+ 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]);
@@ -655,32 +677,41 @@ int main(int argc, char *argv[])
exit(1);
}
}
- for(irow=0;irow<inputReader.nrOfRow()&&!isDifferent;++irow){
+ double rmse=0;
+ // for(irow=0;irow<inputReader.nrOfRow()&&!isDifferent;++irow){
+ for(irow=0;irow<inputReader.nrOfRow();++irow){
//read line in lineInput, lineReference and lineMask
- inputReader.readData(lineInput,GDT_Int16,irow,band_opt[0]);
- if(mask_opt.size())
- maskReader.readData(lineMask,GDT_Int16,irow);
+ inputReader.readData(lineInput,GDT_Float64,irow,band_opt[0]);
double x,y;//geo coordinates
double ireference,jreference;//image coordinates in reference image
+ double imask,jmask;//image coordinates in mask image
for(icol=0;icol<inputReader.nrOfCol();++icol){
//find col in reference
inputReader.image2geo(icol,irow,x,y);
referenceReaderGdal.geo2image(x,y,ireference,jreference);
if(ireference<0||ireference>=referenceReaderGdal.nrOfCol()){
- cerr << ireference << " out of reference range!" << endl;
- cerr << x << " " << y << " " << icol << " " << irow << endl;
- cerr << x << " " << y << " " << ireference << " " << jreference << endl;
- exit(1);
+ if(rmse_opt[0]||regression_opt[0])
+ continue;
+ else{
+ cerr << ireference << " out of reference range!" << endl;
+ cerr << x << " " << y << " " << icol << " " << irow << endl;
+ cerr << x << " " << y << " " << ireference << " " << jreference << endl;
+ exit(1);
+ }
}
if(jreference!=oldreferencerow){
if(jreference<0||jreference>=referenceReaderGdal.nrOfRow()){
- cerr << jreference << " out of reference range!" << endl;
- cerr << x << " " << y << " " << icol << " " << irow << endl;
- cerr << x << " " << y << " " << ireference << " " << jreference << endl;
- exit(1);
+ if(rmse_opt[0]||regression_opt[0])
+ continue;
+ else{
+ cerr << jreference << " out of reference range!" << endl;
+ cerr << x << " " << y << " " << icol << " " << irow << endl;
+ cerr << x << " " << y << " " << ireference << " " << jreference << endl;
+ exit(1);
+ }
}
else{
- referenceReaderGdal.readData(lineReference,GDT_Int16,static_cast<int>(jreference),band_opt[0]);
+ referenceReaderGdal.readData(lineReference,GDT_Float64,static_cast<int>(jreference),band_opt[1]);
oldreferencerow=jreference;
}
}
@@ -694,14 +725,27 @@ int main(int argc, char *argv[])
}
}
if(mask_opt.size()){
- for(int ivalue=0;ivalue<masknodata_opt.size();++ivalue){
- if(lineMask[icol]==masknodata_opt[ivalue]){
- flagged=true;
- break;
- }
+ maskReader.geo2image(x,y,imask,jmask);
+ if(jmask>=0&&jmask<maskReader.nrOfRow()){
+ if(jmask!=oldmaskrow)
+ maskReader.readData(lineMask,GDT_Float64,jmask);
+ for(int ivalue=0;ivalue<msknodata_opt.size();++ivalue){
+ if(lineMask[icol]==msknodata_opt[ivalue]){
+ flagged=true;
+ break;
+ }
+ }
}
}
if(!flagged){
+ if(rmse_opt[0]){//divide by image size to prevent overflow. At the end we need to take care about flagged pixels by normalizing...
+ rmse+=static_cast<double>(lineInput[icol]-lineReference[ireference])*(lineInput[icol]-lineReference[ireference])/inputReader.nrOfCol()/inputReader.nrOfRow();
+ }
+ else if(regression_opt[0]){
+ bufferInput.push_back(lineInput[icol]);
+ bufferReference.push_back(lineReference[ireference]);
+ }
+
if(confusion_opt[0]){
++ntotalValidation;
int rc=distance(referenceRange.begin(),find(referenceRange.begin(),referenceRange.end(),lineReference[ireference]));
@@ -724,22 +768,11 @@ int main(int argc, char *argv[])
}
}
else{//error
- if(output_opt.empty()&&!confusion_opt[0]){
+ if(output_opt.empty()&&!confusion_opt[0]&&!rmse_opt[0]&&!regression_opt[0]){
isDifferent=true;
break;
}
if(output_opt.size()){
- // if(lineInput[icol]<20){//forest
- // if(lineReference[icol]>=20)//gain
- // lineOutput[icol]=lineInput[icol]*10+1;//GAIN is 111,121,131
- // else//forest type changed: mixed
- // lineOutput[icol]=130;//MIXED FOREST
- // }
- // else if(lineReference[icol]<20){//loss
- // lineOutput[icol]=20*10+lineReference[icol];//LOSS is 211 212 213
- // }
- // else//no forest
- // lineOutput[icol]=20*10;//NON FOREST is 200
if(lineInput[icol]>lineReference[ireference])
lineOutput[icol]=valueO_opt[0];//omission error
else
@@ -759,7 +792,7 @@ int main(int argc, char *argv[])
}
if(output_opt.size()){
try{
- gdalWriter.writeData(lineOutput,GDT_Int16,irow);
+ gdalWriter.writeData(lineOutput,GDT_Float64,irow);
}
catch(string errorstring){
cerr << "lineOutput.size(): " << lineOutput.size() << endl;
@@ -768,7 +801,7 @@ int main(int argc, char *argv[])
exit(1);
}
}
- else if(isDifferent&&!confusion_opt[0]){//we can break off here, files are different...
+ else if(isDifferent&&!confusion_opt[0]&&!rmse_opt[0]&&!regression_opt[0]){//we can break off here, files are different...
if(!verbose_opt[0])
pfnProgress(1.0,pszMessage,pProgressArg);
break;
@@ -780,7 +813,38 @@ int main(int argc, char *argv[])
if(output_opt.size())
gdalWriter.close();
else if(!confusion_opt[0]){
- if(isDifferent)
+ if(rmse_opt[0]){
+ double normalization=1.0*inputReader.nrOfCol()*inputReader.nrOfRow()/(inputReader.nrOfCol()*inputReader.nrOfRow()-nflagged);
+ if(verbose_opt[0]){
+ cout << "normalization: " << normalization << endl;
+ cout << "rmse before sqrt and normalization: " << rmse << endl;
+ }
+ cout << "--rmse " << sqrt(rmse/normalization) << endl;
+ }
+ else if(regression_opt[0]){
+ double err=0;
+ double c0=0;
+ double c1=1;
+ statfactory::StatFactory stat;
+ if(bufferInput.size()&&bufferReference.size()){
+ err=stat.linear_regression_err(bufferInput,bufferReference,c0,c1);
+ }
+ if(verbose_opt[0]){
+ cout << "bufferInput.size(): " << bufferInput.size() << endl;
+ cout << "bufferReference.size(): " << bufferReference.size() << endl;
+ double theMin=0;
+ double theMax=0;
+ stat.minmax(bufferInput,bufferInput.begin(),bufferInput.end(),theMin,theMax);
+ cout << "min, max input: " << theMin << ", " << theMax << endl;
+ theMin=0;
+ theMax=0;
+ stat.minmax(bufferReference,bufferReference.begin(),bufferReference.end(),theMin,theMax);
+ cout << "min, max reference: " << theMin << ", " << theMax << endl;
+ }
+ cout << "--c0 " << c0 << "--c1 " << c1 << " --rmse: " << err << endl;
+
+ }
+ else if(isDifferent)
cout << input_opt[0] << " and " << reference_opt[0] << " are different" << endl;
else
cout << input_opt[0] << " and " << reference_opt[0] << " are identical" << endl;
@@ -789,67 +853,9 @@ int main(int argc, char *argv[])
inputReader.close();
if(mask_opt.size())
maskReader.close();
- }
+ }//raster dataset
if(confusion_opt[0]){
- // double totalResult=0;
- // cout << " ";
- // for(int ic=0;ic<inputRange.size();++ic)
- // cout << inputRange[ic] << " ";
- // cout << endl;
- // unsigned int ntotal=0;
- // vector<unsigned int> ntotalclass(referenceRange.size());
- // for(int rc=0;rc<referenceRange.size();++rc){
- // ntotalclass[rc]=0;
- // cout << referenceRange[rc] << " ";
- // //initialize
- // for(int ic=0;ic<inputRange.size();++ic){
- // unsigned int result=0;
- // user[ic]=0;
- // producer[ic]=0;
- // for(int k=0;k<nclass;++k){
- // user[ic]+=resultClass[k][ic];
- // producer[ic]+=resultClass[ic][k];
- // }
- // result=resultClass[rc][ic];
- // ntotal+=result;
- // ntotalclass[rc]+=result;
- // if(ic==rc){
- // totalResult+=result;
- // }
- // cout << result << " ";
- // }
- // cout << endl;
- // }
- // if(verbose_opt[0]){
- // cout << "totalResult: " << totalResult << endl;
- // cout << "ntotalValidation: " << ntotalValidation << endl;
- // cout << "nflagged: " << nflagged << endl;
- // cout << "ntotal: " << ntotal << endl;
- // }
- // totalResult*=100.0/ntotal;
- // int nclass0=0;//number of classes without any reference
- // for(int rc=0;rc<referenceRange.size();++rc)
- // if(!nvalidation[rc])
- // ++nclass0;
- // double pChance=0;
- // double pCorrect=0;
- // double totalEntries=0;
- // cout << "class #samples userAcc prodAcc" << endl;
- // for(int rc=0;rc<referenceRange.size();++rc){
- // totalEntries+=user[rc];
- // pChance+=user[rc]*producer[rc];
- // pCorrect+=resultClass[rc][rc];
- // cout << referenceRange[rc] << " " << nvalidation[rc] << " ";
- // cout << 100.0*resultClass[rc][rc]/user[rc] << " "//user accuracy
- // << 100.0*resultClass[rc][rc]/producer[rc] << " " << endl;//producer accuracy
- // }
- // pCorrect/=totalEntries;
- // pChance/=totalEntries*totalEntries;
- // double kappa=pCorrect-pChance;
- // kappa/=1-pChance;
- // cout << "Kappa: " << kappa << endl;
- // cout << "total weighted: " << static_cast<int>(0.5+totalResult) << endl;
assert(cm.nReference());
cout << cm << endl;
cout << "class #samples userAcc prodAcc" << endl;
diff --git a/src/apps/pkdsm2shadow.cc b/src/apps/pkdsm2shadow.cc
index ea8d586..b91a900 100644
--- a/src/apps/pkdsm2shadow.cc
+++ b/src/apps/pkdsm2shadow.cc
@@ -41,13 +41,16 @@ int main(int argc,char **argv) {
Optionpk<double> sza_opt("sza", "sza", "Sun zenith angle.");
Optionpk<double> saa_opt("saa", "saa", "Sun azimuth angle (N=0 E=90 S=180 W=270).");
Optionpk<int> flag_opt("f", "flag", "Flag to put in image if pixel shadow", 0);
- Optionpk<double> scale_opt("s", "scale", "scale used for input dsm: height=scale*input+offset");
- Optionpk<double> offset_opt("off", "offset", "offset used for input dsm: height=scale*input+offset");
Optionpk<std::string> otype_opt("ot", "otype", "Data type for output image ({Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/CInt16/CInt32/CFloat32/CFloat64}). Empty string: inherit type from input image", "");
Optionpk<string> oformat_opt("of", "oformat", "Output image format (see also gdal_translate). Empty string: inherit from input image");
Optionpk<string> colorTable_opt("ct", "ct", "color table (file with 5 columns: id R G B ALFA (0: transparent, 255: solid)");
Optionpk<string> option_opt("co", "co", "Creation option for output file. Multiple options can be specified.");
- Optionpk<short> verbose_opt("v", "verbose", "verbose mode if > 0", 0);
+ Optionpk<double> scale_opt("s", "scale", "scale used for input dsm: height=scale*input+offset");
+ Optionpk<double> offset_opt("off", "offset", "offset used for input dsm: height=scale*input+offset");
+ Optionpk<short> verbose_opt("v", "verbose", "verbose mode if > 0", 0,2);
+
+ scale_opt.setHide(1);
+ offset_opt.setHide(1);
bool doProcess;//stop process when program was invoked with help option (-h --help)
try{
@@ -69,6 +72,9 @@ int main(int argc,char **argv) {
exit(0);
}
if(!doProcess){
+ cout << endl;
+ cout << "Usage: pkdsm2shadow -i input.txt -o output [-sza angle] [-saa angle]" << 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
}
diff --git a/src/apps/pkdumpimg.cc b/src/apps/pkdumpimg.cc
index e763d86..c6334e6 100644
--- a/src/apps/pkdumpimg.cc
+++ b/src/apps/pkdumpimg.cc
@@ -49,8 +49,14 @@ int main(int argc, char *argv[])
Optionpk<string> resample_opt("r", "resampling-method", "Resampling method (near: nearest neighbour, bilinear: bi-linear interpolation).", "near");
Optionpk<short> dstnodata_opt("dstnodata", "dstnodata", "nodata value for output if out of bounds.", 0);
Optionpk<double> srcnodata_opt("srcnodata", "srcnodata", "set no data value(s) for input image");
- Optionpk<short> verbose_opt("v", "verbose", "verbose (Default: 0)", 0);
+ Optionpk<short> verbose_opt("v", "verbose", "verbose (Default: 0)", 0,2);
+ dx_opt.setHide(1);
+ dy_opt.setHide(1);
+ resample_opt.setHide(1);
+ srcnodata_opt.setHide(1);
+ dstnodata_opt.setHide(1);
+
bool doProcess;//stop process when program was invoked with help option (-h --help)
try{
doProcess=input_opt.retrieveOption(argc,argv);
@@ -74,6 +80,9 @@ int main(int argc, char *argv[])
exit(0);
}
if(!doProcess){
+ cout << endl;
+ cout << "Usage: pkdumpimg -i input.txt [-o output]" << 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
}
diff --git a/src/apps/pkdumpogr.cc b/src/apps/pkdumpogr.cc
index 2a4e305..ccdbccf 100644
--- a/src/apps/pkdumpogr.cc
+++ b/src/apps/pkdumpogr.cc
@@ -33,7 +33,7 @@ int main(int argc, char *argv[])
Optionpk<string> attribute_opt("n", "name", "names of the attributes to select. Each attribute is stored in a separate band. Default is ALL: write all attributes", "ALL");
Optionpk<bool> pos_opt("pos","pos","include position (x and y)",false);
Optionpk<bool> transpose_opt("t","transpose","transpose output (does not work for -n ALL ",false);
- Optionpk<short> verbose_opt("v", "verbose", "verbose (Default: 0)", 0);
+ Optionpk<short> verbose_opt("v", "verbose", "verbose (Default: 0)", 0,2);
bool doProcess;//stop process when program was invoked with help option (-h --help)
try{
@@ -50,6 +50,9 @@ int main(int argc, char *argv[])
exit(0);
}
if(!doProcess){
+ cout << endl;
+ cout << "Usage: pkdumpogr -i input.txt [-o output]" << 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
}
diff --git a/src/apps/pkextract.cc b/src/apps/pkextract.cc
index e65124f..6aea515 100644
--- a/src/apps/pkextract.cc
+++ b/src/apps/pkextract.cc
@@ -54,7 +54,7 @@ int main(int argc, char *argv[])
Optionpk<string> ftype_opt("ft", "ftype", "Field type (only Real or Integer)", "Real");
Optionpk<string> ltype_opt("lt", "ltype", "Label type: In16 or String", "Integer");
Optionpk<bool> polygon_opt("polygon", "polygon", "Create OGRPolygon as geometry instead of OGRPoint.", false);
- Optionpk<int> band_opt("b", "band", "Band index(es) to extract. Leave empty to use all bands");
+ Optionpk<int> band_opt("b", "band", "Band index(es) to extract (0 based). Leave empty to use all bands");
Optionpk<string> rule_opt("r", "rule", "Rule how to report image information per feature (only for vector sample). point (value at each point or at centroid if polygon), centroid, mean, stdev, median, proportion, min, max, mode, sum.", "centroid");
Optionpk<double> srcnodata_opt("srcnodata", "srcnodata", "Invalid value(s) for input image");
Optionpk<int> bndnodata_opt("bndnodata", "bndnodata", "Band(s) in input image to check if pixel is valid (used for srcnodata)", 0);
@@ -66,7 +66,18 @@ int main(int argc, char *argv[])
Optionpk<short> down_opt("down", "down", "Down sampling factor (for raster sample datasets only). Can be used to create grid points", 1);
Optionpk<short> buffer_opt("buf", "buffer", "Buffer for calculating statistics for point features ");
Optionpk<bool> disc_opt("circ", "circular", "Use a circular disc kernel buffer (for vector point sample datasets only, use in combination with buffer option)", false);
- Optionpk<short> verbose_opt("v", "verbose", "Verbose mode if > 0", 0);
+ Optionpk<short> verbose_opt("v", "verbose", "Verbose mode if > 0", 0,2);
+
+ bndnodata_opt.setHide(1);
+ srcnodata_opt.setHide(1);
+ polythreshold_opt.setHide(1);
+ test_opt.setHide(1);
+ fieldname_opt.setHide(1);
+ label_opt.setHide(1);
+ geo_opt.setHide(1);
+ down_opt.setHide(1);
+ buffer_opt.setHide(1);
+ disc_opt.setHide(1);
bool doProcess;//stop process when program was invoked with help option (-h --help)
try{
@@ -103,6 +114,9 @@ int main(int argc, char *argv[])
exit(0);
}
if(!doProcess){
+ cout << endl;
+ cout << "Usage: pkextract -i input [-s sample | -rand number | -grid size] -o output" << 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
}
@@ -131,6 +145,7 @@ int main(int argc, char *argv[])
if(verbose_opt[0])
std::cout << class_opt << std::endl;
statfactory::StatFactory stat;
+ stat.setNoDataValues(srcnodata_opt);
Vector2d<unsigned int> posdata;
unsigned long int nsample=0;
unsigned long int ntotalvalid=0;
@@ -392,12 +407,16 @@ int main(int argc, char *argv[])
imgReader.readData(imgBuffer[iband],GDT_Float64,static_cast<int>(jimg),theBand);
assert(imgBuffer[iband].size()==imgReader.nrOfCol());
if(srcnodata_opt.size()){
- vector<int>::const_iterator bndit=find(bndnodata_opt.begin(),bndnodata_opt.end(),theBand);
- if(bndit!=bndnodata_opt.end()){
- vector<int>::const_iterator bndit=find(bndnodata_opt.begin(),bndnodata_opt.end(),theBand);
- if(bndit!=bndnodata_opt.end()){
- if(imgBuffer[iband][static_cast<int>(iimg)]==srcnodata_opt[theBand])
- valid=false;
+ vector<int>::const_iterator bndit=bndnodata_opt.begin();
+ vector<double>::const_iterator srcit=srcnodata_opt.begin();
+ while(bndit!=bndnodata_opt.end()&&srcit!=srcnodata_opt.end()){
+ if((*bndit==theBand)&&(*srcit==imgBuffer[iband][static_cast<int>(iimg)])){
+ valid=false;
+ break;
+ }
+ else{
+ ++bndit;
+ ++srcit;
}
}
}
@@ -588,11 +607,19 @@ int main(int argc, char *argv[])
int theBand=(band_opt.size()) ? band_opt[iband] : iband;
imgReader.readData(imgBuffer[iband],GDT_Float64,static_cast<int>(jimg),theBand);
assert(imgBuffer[iband].size()==imgReader.nrOfCol());
+
if(srcnodata_opt.size()){
- vector<int>::const_iterator bndit=find(bndnodata_opt.begin(),bndnodata_opt.end(),theBand);
- if(bndit!=bndnodata_opt.end()){
- if(imgBuffer[iband][static_cast<int>(iimg)]==srcnodata_opt[theBand])
+ vector<int>::const_iterator bndit=bndnodata_opt.begin();
+ vector<double>::const_iterator srcit=srcnodata_opt.begin();
+ while(bndit!=bndnodata_opt.end()&&srcit!=srcnodata_opt.end()){
+ if((*bndit==theBand)&&(*srcit==imgBuffer[iband][static_cast<int>(iimg)])){
valid=false;
+ break;
+ }
+ else{
+ ++bndit;
+ ++srcit;
+ }
}
}
}
@@ -1536,13 +1563,13 @@ int main(int argc, char *argv[])
writePointFeature = OGRFeature::CreateFeature(writeLayer->GetLayerDefn());
if(verbose_opt[0]>1)
std::cout << "copying fields from polygons " << std::endl;
- //test
- // writePointFeature->SetGeometry(&thePoint);
if(writePointFeature->SetFrom(readFeature)!= OGRERR_NONE)
cerr << "writing feature failed" << std::endl;
+ if(verbose_opt[0]>1)
+ std::cout << "set geometry as point " << std::endl;
//test
writePointFeature->SetGeometry(&thePoint);
- assert(wkbFlatten(writePolygonFeature->GetGeometryRef()->getGeometryType()) == wkbPoint);
+ assert(wkbFlatten(writePointFeature->GetGeometryRef()->getGeometryType()) == wkbPoint);
//test
// OGRGeometry *updateGeometry;
// updateGeometry = writePointFeature->GetGeometryRef();
@@ -2027,7 +2054,7 @@ int main(int argc, char *argv[])
cerr << "writing feature failed" << std::endl;
//test
writePointFeature->SetGeometry(&thePoint);
- assert(wkbFlatten(writePolygonFeature->GetGeometryRef()->getGeometryType()) == wkbPoint);
+ assert(wkbFlatten(writePointFeature->GetGeometryRef()->getGeometryType()) == wkbPoint);
// OGRGeometry *updateGeometry;
// updateGeometry = writePointFeature->GetGeometryRef();
// OGRPoint *poPoint = (OGRPoint *) updateGeometry;
diff --git a/src/apps/pkfillnodata.cc b/src/apps/pkfillnodata.cc
index 1f003a8..681d1bf 100644
--- a/src/apps/pkfillnodata.cc
+++ b/src/apps/pkfillnodata.cc
@@ -26,6 +26,8 @@ extern "C" {
#include <string>
#include "base/Optionpk.h"
+using namespace std;
+
int main(int argc,char **argv) {
Optionpk<std::string> input_opt("i", "input", "Input raster dataset");
Optionpk<int> band_opt("b", "band", "band(s) to process (Default is -1: process all bands)");
@@ -33,7 +35,10 @@ int main(int argc,char **argv) {
Optionpk<std::string> output_opt("o", "output", "Output image file");
Optionpk<double> distance_opt("d", "distance", "Maximum number of pixels to search in all directions to find values to interpolate from", 0);
Optionpk<int> iteration_opt("it", "iteration", "Number of 3x3 smoothing filter passes to run (default 0)", 0);
- Optionpk<short> verbose_opt("v", "verbose", "verbose", 0);
+ Optionpk<short> verbose_opt("v", "verbose", "verbose", 0,2);
+
+ distance_opt.setHide(1);
+ iteration_opt.setHide(1);
bool doProcess;//stop process when program was invoked with help option (-h --help)
try{
@@ -50,6 +55,9 @@ int main(int argc,char **argv) {
exit(0);
}
if(!doProcess){
+ cout << endl;
+ cout << "Usage: pkfillnodata -i input.txt -m mask -o output" << 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
}
diff --git a/src/apps/pkfilter.cc b/src/apps/pkfilter.cc
index 67489ad..3c352ce 100644
--- a/src/apps/pkfilter.cc
+++ b/src/apps/pkfilter.cc
@@ -39,27 +39,31 @@ using namespace std;
int main(int argc,char **argv) {
Optionpk<std::string> input_opt("i","input","input image file");
Optionpk<std::string> output_opt("o", "output", "Output image file");
- Optionpk<std::string> tmpdir_opt("tmp", "tmp", "Temporary directory","/tmp",2);
+ // Optionpk<std::string> tmpdir_opt("tmp", "tmp", "Temporary directory","/tmp",2);
Optionpk<bool> disc_opt("circ", "circular", "circular disc kernel for dilation and erosion", false);
// Optionpk<double> angle_opt("a", "angle", "angle used for directional filtering in dilation (North=0, East=90, South=180, West=270).");
- Optionpk<std::string> method_opt("f", "filter", "filter function (median, var, min, max, sum, mean, dilate, erode, close, open, homog (central pixel must be identical to all other pixels within window), heterog, sobelx (horizontal edge detection), sobely (vertical edge detection), sobelxy (diagonal edge detection NE-SW),sobelyx (diagonal edge detection NW-SE), smooth, density, countid, mode (majority voting, only for classes), smoothnodata (smooth nodata values only) values, threshold [...]
+ Optionpk<std::string> method_opt("f", "filter", "filter function (median, var, min, max, sum, mean, dilate, erode, close, open, homog (central pixel must be identical to all other pixels within window), heterog, sobelx (horizontal edge detection), sobely (vertical edge detection), sobelxy (diagonal edge detection NE-SW),sobelyx (diagonal edge detection NW-SE), smooth, density, countid, mode (majority voting, only for classes), smoothnodata (smooth nodata values only) values, threshold [...]
Optionpk<std::string> resample_opt("r", "resampling-method", "Resampling method for shifting operation (near: nearest neighbour, bilinear: bi-linear interpolation).", "near");
Optionpk<double> dimX_opt("dx", "dx", "filter kernel size in x, better use odd value to avoid image shift", 3);
Optionpk<double> dimY_opt("dy", "dy", "filter kernel size in y, better use odd value to avoid image shift", 3);
Optionpk<int> dimZ_opt("dz", "dz", "filter kernel size in z (band or spectral dimension), must be odd (example: 3).. Set dz>0 if 1-D filter must be used in band domain");
Optionpk<std::string> wavelet_type_opt("wt", "wavelet", "wavelet type: daubechies,daubechies_centered, haar, haar_centered, bspline, bspline_centered", "daubechies");
Optionpk<int> family_opt("wf", "family", "wavelet family (vanishing moment, see also http://www.gnu.org/software/gsl/manual/html_node/DWT-Initialization.html)", 4);
+ Optionpk<int> savgolay_nl_opt("nl", "nl", "Number of leftward (past) data points used in Savitzky-Golay filter)", 2);
+ Optionpk<int> savgolay_nr_opt("nr", "nr", "Number of rightward (future) data points used in Savitzky-Golay filter)", 2);
+ Optionpk<int> savgolay_ld_opt("ld", "ld", "order of the derivative desired in Savitzky-Golay filter (e.g., ld=0 for smoothed function)", 0);
+ Optionpk<int> savgolay_m_opt("m", "m", "order of the smoothing polynomial in Savitzky-Golay filter, also equal to the highest conserved moment; usual values are m = 2 or m = 4)", 2);
Optionpk<short> class_opt("class", "class", "class value(s) to use for density, erosion, dilation, openening and closing, thresholding");
Optionpk<double> threshold_opt("t", "threshold", "threshold value(s) to use for threshold filter (one for each class), or threshold to cut for dwt_cut (use 0 to keep all) or dwt_cut_from, or sigma for shift", 0);
Optionpk<double> nodata_opt("nodata", "nodata", "nodata value(s) (used for smoothnodata filter)");
Optionpk<std::string> tap_opt("tap", "tap", "text file containing taps used for spatial filtering (from ul to lr). Use dimX and dimY to specify tap dimensions in x and y. Leave empty for not using taps");
Optionpk<double> tapz_opt("tapz", "tapz", "taps used for spectral filtering");
- Optionpk<string> padding_opt("pad","pad", "Padding method for filtering (how to handle edge effects). Choose between: symmetric, replicate, circular, constant (pad with 0).", "symmetric");
+ Optionpk<string> padding_opt("pad","pad", "Padding method for filtering (how to handle edge effects). Choose between: symmetric, replicate, circular, zero (pad with 0).", "symmetric");
Optionpk<double> fwhm_opt("fwhm", "fwhm", "list of full width half to apply spectral filtering (-fwhm band1 -fwhm band2 ...)");
Optionpk<std::string> srf_opt("srf", "srf", "list of ASCII files containing spectral response functions (two columns: wavelength response)");
Optionpk<double> wavelengthIn_opt("win", "wavelengthIn", "list of wavelengths in input spectrum (-win band1 -win band2 ...)");
Optionpk<double> wavelengthOut_opt("wout", "wavelengthOut", "list of wavelengths in output spectrum (-wout band1 -wout band2 ...)");
- Optionpk<std::string> interpolationType_opt("interp", "interp", "type of interpolation for spectral filtering (see http://www.gnu.org/software/gsl/manual/html_node/Interpolation-Types.html)","akima",1);
+ Optionpk<std::string> interpolationType_opt("interp", "interp", "type of interpolation for spectral filtering (see http://www.gnu.org/software/gsl/manual/html_node/Interpolation-Types.html)","akima");
Optionpk<std::string> otype_opt("ot", "otype", "Data type for output image ({Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/CInt16/CInt32/CFloat32/CFloat64}). Empty string: inherit type from input image","");
Optionpk<string> oformat_opt("of", "oformat", "Output image format (see also gdal_translate). Empty string: inherit from input image");
Optionpk<string> colorTable_opt("ct", "ct", "color table (file with 5 columns: id R G B ALFA (0: transparent, 255: solid). Use none to ommit color table");
@@ -71,31 +75,60 @@ int main(int argc,char **argv) {
Optionpk<bool> l2_opt("l2","l2", "obtain shortest object length for linear feature",false,2);
Optionpk<bool> a1_opt("a1","a1", "obtain angle found for longest object length for linear feature",false);
Optionpk<bool> a2_opt("a2","a2", "obtain angle found for shortest object length for linear feature",false);
- Optionpk<short> verbose_opt("v", "verbose", "verbose mode if > 0", 0);
+ Optionpk<short> verbose_opt("v", "verbose", "verbose mode if > 0", 0,2);
+
+ resample_opt.setHide(1);
+ option_opt.setHide(1);
+ wavelet_type_opt.setHide(1);
+ savgolay_nl_opt.setHide(1);
+ savgolay_nr_opt.setHide(1);
+ savgolay_ld_opt.setHide(1);
+ savgolay_m_opt.setHide(1);
+ class_opt.setHide(1);
+ threshold_opt.setHide(1);
+ tap_opt.setHide(1);
+ tapz_opt.setHide(1);
+ padding_opt.setHide(1);
+ wavelengthIn_opt.setHide(1);
+ wavelengthOut_opt.setHide(1);
+ down_opt.setHide(1);
+ beta_opt.setHide(1);
+ eps_opt.setHide(1);
+ l1_opt.setHide(1);
+ l2_opt.setHide(1);
+ a1_opt.setHide(1);
+ a2_opt.setHide(1);
+ interpolationType_opt.setHide(1);
+ otype_opt.setHide(1);
+ oformat_opt.setHide(1);
+ colorTable_opt.setHide(1);
+ disc_opt.setHide(1);
bool doProcess;//stop process when program was invoked with help option (-h --help)
try{
doProcess=input_opt.retrieveOption(argc,argv);
output_opt.retrieveOption(argc,argv);
- tmpdir_opt.retrieveOption(argc,argv);
- disc_opt.retrieveOption(argc,argv);
+ // tmpdir_opt.retrieveOption(argc,argv);
// angle_opt.retrieveOption(argc,argv);
method_opt.retrieveOption(argc,argv);
- resample_opt.retrieveOption(argc,argv);
+ srf_opt.retrieveOption(argc,argv);
+ fwhm_opt.retrieveOption(argc,argv);
dimX_opt.retrieveOption(argc,argv);
dimY_opt.retrieveOption(argc,argv);
dimZ_opt.retrieveOption(argc,argv);
+ nodata_opt.retrieveOption(argc,argv);
+ resample_opt.retrieveOption(argc,argv);
option_opt.retrieveOption(argc,argv);
wavelet_type_opt.retrieveOption(argc,argv);
- family_opt.retrieveOption(argc,argv);
+ savgolay_nl_opt.retrieveOption(argc,argv);
+ savgolay_nr_opt.retrieveOption(argc,argv);
+ savgolay_ld_opt.retrieveOption(argc,argv);
+ savgolay_m_opt.retrieveOption(argc,argv);
class_opt.retrieveOption(argc,argv);
threshold_opt.retrieveOption(argc,argv);
- nodata_opt.retrieveOption(argc,argv);
tap_opt.retrieveOption(argc,argv);
tapz_opt.retrieveOption(argc,argv);
padding_opt.retrieveOption(argc,argv);
- fwhm_opt.retrieveOption(argc,argv);
- srf_opt.retrieveOption(argc,argv);
wavelengthIn_opt.retrieveOption(argc,argv);
wavelengthOut_opt.retrieveOption(argc,argv);
down_opt.retrieveOption(argc,argv);
@@ -109,6 +142,7 @@ int main(int argc,char **argv) {
otype_opt.retrieveOption(argc,argv);
oformat_opt.retrieveOption(argc,argv);
colorTable_opt.retrieveOption(argc,argv);
+ disc_opt.retrieveOption(argc,argv);
verbose_opt.retrieveOption(argc,argv);
}
catch(string predefinedString){
@@ -116,6 +150,9 @@ int main(int argc,char **argv) {
exit(0);
}
if(!doProcess){
+ cout << endl;
+ cout << "Usage: pkfilter -i input -o ouptut [-f filter | -srf file [-srf file]* | -fwhm value [-fwhm value]*]" << 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
}
@@ -192,9 +229,19 @@ int main(int argc,char **argv) {
}
else{
int nband=input.nrOfBand();
- if(dimZ_opt.size())
- if(dimZ_opt[0]==1)
+ if(dimZ_opt.size()){
+ if(filter2d::Filter2d::getFilterType(method_opt[0])==filter2d::smoothnodata)
+ dimZ_opt[0]=nband;
+ if(dimZ_opt[0]==1){
nband=1;
+ if(verbose_opt[0])
+ std::cout << "opening single band output image " << output_opt[0] << std::endl;
+ }
+ else if(verbose_opt[0])
+ std::cout << "opening multi-band output image " << output_opt[0] << std::endl;
+ }
+ else
+ std::cout << "opening output image " << output_opt[0] << std::endl;
output.open(output_opt[0],(input.nrOfCol()+down_opt[0]-1)/down_opt[0],(input.nrOfRow()+down_opt[0]-1)/down_opt[0],nband,theType,imageType,option_opt);
}
}
@@ -287,7 +334,7 @@ int main(int argc,char **argv) {
std::cout<< tapz_opt[itap] << " ";
std::cout<< std::endl;
}
- filter1d.setTaps(tapz_opt);
+ filter1d.setTaps(tapz_opt);
filter1d.filter(input,output);
}
else if(fwhm_opt.size()){
@@ -420,10 +467,12 @@ int main(int argc,char **argv) {
std::cerr << "Error: down option not supported for morphological operator" << std::endl;
exit(1);
}
- ostringstream tmps;
- tmps << tmpdir_opt[0] << "/dilation_" << getpid() << ".tif";
+
+ // ostringstream tmps;
+ // tmps << tmpdir_opt[0] << "/dilation_" << getpid() << ".tif";
ImgWriterGdal tmpout;
- tmpout.open(tmps.str(),input);
+ tmpout.open("/vsimem/dilation",input);
+ // tmpout.open(tmps.str(),input);
try{
if(dimZ_opt.size()){
filter1d.morphology(input,tmpout,"dilate",dimZ_opt[0]);
@@ -436,19 +485,21 @@ int main(int argc,char **argv) {
std::cout<< errorString;
exit(1);
}
- tmpout.close();
+ // tmpout.close();
ImgReaderGdal tmpin;
- tmpin.open(tmps.str());
+ tmpin.open("/vsimem/dilation");
+ // tmpin.open(tmps.str());
if(dimZ_opt.size()){
filter1d.morphology(tmpin,output,"erode",dimZ_opt[0]);
}
else{
- filter2d.morphology(tmpin,output,"erode",dimX_opt[0],dimY_opt[0],angle_opt,disc_opt[0]);
+ filter2d.morphology(tmpin,output,"erode",dimX_opt[0],dimY_opt[0],angle_opt,disc_opt[0]);
}
tmpin.close();
- if(remove(tmps.str().c_str( )) !=0){
- cerr << "could not remove " << tmps.str() << std::endl;
- }
+ tmpout.close();
+ // if(remove(tmps.str().c_str( )) !=0){
+ // cerr << "could not remove " << tmps.str() << std::endl;
+ // }
break;
}
case(filter2d::open):{//opening
@@ -456,19 +507,27 @@ int main(int argc,char **argv) {
std::cerr << "Error: down option not supported for morphological operator" << std::endl;
exit(1);
}
- ostringstream tmps;
- tmps << tmpdir_opt[0] << "/erosion_" << getpid() << ".tif";
+ // ostringstream tmps;
+ // tmps << tmpdir_opt[0] << "/erosion_" << getpid() << ".tif";
ImgWriterGdal tmpout;
- tmpout.open(tmps.str(),input);
- if(dimZ_opt.size()){
- filter1d.morphology(input,tmpout,"erode",dimZ_opt[0]);
+ tmpout.open("/vsimem/erosion",input);
+ // tmpout.open(tmps.str(),input);
+ try{
+ if(dimZ_opt.size()){
+ filter1d.morphology(input,tmpout,"erode",dimZ_opt[0]);
+ }
+ else{
+ filter2d.morphology(input,tmpout,"erode",dimX_opt[0],dimY_opt[0],angle_opt,disc_opt[0]);
+ }
}
- else{
- filter2d.morphology(input,tmpout,"erode",dimX_opt[0],dimY_opt[0],angle_opt,disc_opt[0]);
+ catch(std::string errorString){
+ std::cout<< errorString;
+ exit(1);
}
- tmpout.close();
+ // tmpout.close();
ImgReaderGdal tmpin;
- tmpin.open(tmps.str());
+ tmpin.open("/vsimem/erosion");
+ // tmpin.open(tmps.str());
if(dimZ_opt.size()){
filter1d.morphology(tmpin,output,"dilate",dimZ_opt[0]);
}
@@ -476,9 +535,10 @@ int main(int argc,char **argv) {
filter2d.morphology(tmpin,output,"dilate",dimX_opt[0],dimY_opt[0],angle_opt,disc_opt[0]);
}
tmpin.close();
- if(remove(tmps.str().c_str( )) !=0){
- cerr << "could not remove " << tmps.str() << std::endl;
- }
+ tmpout.close();
+ // if(remove(tmps.str().c_str( )) !=0){
+ // cerr << "could not remove " << tmps.str() << std::endl;
+ // }
break;
}
case(filter2d::homog):{//spatially homogeneous
@@ -654,7 +714,13 @@ int main(int argc,char **argv) {
std::cerr << "Error: down option not supported for this filter" << std::endl;
exit(1);
}
- filter2d.smoothNoData(input,output,dimX_opt[0],dimY_opt[0]);
+ if(dimZ_opt.size()){
+ if(verbose_opt[0])
+ std::cout<< "1-D filtering: smooth" << std::endl;
+ filter1d.smoothNoData(input,interpolationType_opt[0],output);
+ }
+ else
+ filter2d.smoothNoData(input,output,dimX_opt[0],dimY_opt[0]);
break;
}
case(filter2d::dwt):
@@ -711,6 +777,24 @@ int main(int argc,char **argv) {
exit(1);
}
break;
+ case(filter2d::savgolay):{
+ assert(savgolay_nl_opt.size());
+ assert(savgolay_nr_opt.size());
+ assert(savgolay_ld_opt.size());
+ assert(savgolay_m_opt.size());
+ if(verbose_opt[0])
+ std::cout << "Calculating Savitzky-Golay coefficients: " << endl;
+ filter1d.getSavGolayCoefficients(tapz_opt, input.nrOfBand(), savgolay_nl_opt[0], savgolay_nr_opt[0], savgolay_ld_opt[0], savgolay_m_opt[0]);
+ if(verbose_opt[0]){
+ std::cout << "taps (size is " << tapz_opt.size() << "): ";
+ for(int itap=0;itap<tapz_opt.size();++itap)
+ std::cout<< tapz_opt[itap] << " ";
+ std::cout<< std::endl;
+ }
+ filter1d.setTaps(tapz_opt);
+ filter1d.filter(input,output);
+ break;
+ }
case(filter2d::threshold):
filter2d.setThresholds(threshold_opt);//deliberate fall through
case(filter2d::density):
diff --git a/src/apps/pkfilterascii.cc b/src/apps/pkfilterascii.cc
index 7e43e55..9699244 100644
--- a/src/apps/pkfilterascii.cc
+++ b/src/apps/pkfilterascii.cc
@@ -54,7 +54,18 @@ int main(int argc,char **argv) {
Optionpk<double> wavelengthOut_opt("wout", "wavelengthOut", "list of wavelengths in output spectrum (-wout band1 -wout band2 ...)");
Optionpk<std::string> interpolationType_opt("interp", "interp", "type of interpolation for spectral filtering (see http://www.gnu.org/software/gsl/manual/html_node/Interpolation-Types.html)","akima");
Optionpk<bool> transpose_opt("t", "transpose", "transpose output with samples in rows and wavelengths in cols", false);
- Optionpk<short> verbose_opt("v", "verbose", "verbose mode if > 0", 0);
+ Optionpk<short> verbose_opt("v", "verbose", "verbose mode if > 0", 0,2);
+
+ tapZ_opt.setHide(1);
+ fwhm_opt.setHide(1);
+ srf_opt.setHide(1);
+ wavelengthIn_opt.setHide(1);
+ wavelengthOut_opt.setHide(1);
+ interpolationType_opt.setHide(1);
+ transpose_opt.setHide(1);
+ wavelet_type_opt.setHide(1);
+ family_opt.setHide(1);
+ threshold_opt.setHide(1);
bool doProcess;//stop process when program was invoked with help option (-h --help)
try{
@@ -62,9 +73,6 @@ int main(int argc,char **argv) {
output_opt.retrieveOption(argc,argv);
inputCols_opt.retrieveOption(argc,argv);
method_opt.retrieveOption(argc,argv);
- wavelet_type_opt.retrieveOption(argc,argv);
- family_opt.retrieveOption(argc,argv);
- threshold_opt.retrieveOption(argc,argv);
dimZ_opt.retrieveOption(argc,argv);
tapZ_opt.retrieveOption(argc,argv);
fwhm_opt.retrieveOption(argc,argv);
@@ -73,6 +81,9 @@ int main(int argc,char **argv) {
wavelengthOut_opt.retrieveOption(argc,argv);
interpolationType_opt.retrieveOption(argc,argv);
transpose_opt.retrieveOption(argc,argv);
+ wavelet_type_opt.retrieveOption(argc,argv);
+ family_opt.retrieveOption(argc,argv);
+ threshold_opt.retrieveOption(argc,argv);
verbose_opt.retrieveOption(argc,argv);
}
catch(string predefinedString){
@@ -80,6 +91,9 @@ int main(int argc,char **argv) {
exit(0);
}
if(!doProcess){
+ cout << endl;
+ cout << "Usage: pkfilterascii -i input.txt [-ic 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
}
diff --git a/src/apps/pkfilterdem.cc b/src/apps/pkfilterdem.cc
index c2448dd..8848258 100644
--- a/src/apps/pkfilterdem.cc
+++ b/src/apps/pkfilterdem.cc
@@ -33,7 +33,7 @@ using namespace std;
int main(int argc,char **argv) {
Optionpk<std::string> input_opt("i","input","input image file");
Optionpk<std::string> output_opt("o", "output", "Output image file");
- Optionpk<std::string> tmpdir_opt("tmp", "tmp", "Temporary directory","/tmp",2);
+ // Optionpk<std::string> tmpdir_opt("tmp", "tmp", "Temporary directory","/tmp",2);
Optionpk<bool> disc_opt("circ", "circular", "circular disc kernel for dilation and erosion", false);
Optionpk<string> postFilter_opt("f", "filter", "post processing filter: vito, etew_min, promorph (progressive morphological filter),open,close).");
Optionpk<double> dim_opt("dim", "dim", "maximum filter kernel size", 17);
@@ -45,16 +45,25 @@ int main(int argc,char **argv) {
Optionpk<string> colorTable_opt("ct", "ct", "color table (file with 5 columns: id R G B ALFA (0: transparent, 255: solid). Use none to ommit color table");
Optionpk<string> option_opt("co", "co", "Creation option for output file. Multiple options can be specified.");
Optionpk<short> nodata_opt("nodata", "nodata", "nodata value");
- Optionpk<short> verbose_opt("v", "verbose", "verbose mode if > 0", 0);
+ Optionpk<short> verbose_opt("v", "verbose", "verbose mode if > 0", 0,2);
+
+ disc_opt.setHide(1);
+ maxSlope_opt.setHide(1);
+ hThreshold_opt.setHide(1);
+ minChange_opt.setHide(1);
+ otype_opt.setHide(1);
+ oformat_opt.setHide(1);
+ colorTable_opt.setHide(1);
+ nodata_opt.setHide(1);
bool doProcess;//stop process when program was invoked with help option (-h --help)
try{
doProcess=input_opt.retrieveOption(argc,argv);
output_opt.retrieveOption(argc,argv);
- tmpdir_opt.retrieveOption(argc,argv);
- disc_opt.retrieveOption(argc,argv);
+ // tmpdir_opt.retrieveOption(argc,argv);
postFilter_opt.retrieveOption(argc,argv);
dim_opt.retrieveOption(argc,argv);
+ disc_opt.retrieveOption(argc,argv);
maxSlope_opt.retrieveOption(argc,argv);
hThreshold_opt.retrieveOption(argc,argv);
minChange_opt.retrieveOption(argc,argv);
@@ -69,6 +78,9 @@ int main(int argc,char **argv) {
exit(0);
}
if(!doProcess){
+ cout << endl;
+ cout << "Usage: pkfilterdem -i input.txt -o output" << 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
}
diff --git a/src/apps/pkfsann.cc b/src/apps/pkfsann.cc
index 08946e5..3d931af 100644
--- a/src/apps/pkfsann.cc
+++ b/src/apps/pkfsann.cc
@@ -1,5 +1,5 @@
/**********************************************************************
-pkfsann.cc: feature selection for nn classifier
+pkfsann.cc: feature selection for artificial neural network classifier pkann
Copyright (C) 2008-2014 Pieter Kempeneers
This file is part of pktools
@@ -196,7 +196,30 @@ int main(int argc, char *argv[])
Optionpk<float> weights_opt("w", "weights", "weights for neural network. Apply to fully connected network only, starting from first input neuron to last output neuron, including the bias neurons (last neuron in each but last layer)", 0.0);
Optionpk<float> learning_opt("l", "learning", "learning rate (default: 0.7)", 0.7);
Optionpk<unsigned int> maxit_opt("\0", "maxit", "number of maximum iterations (epoch) (default: 500)", 500);
- Optionpk<short> verbose_opt("v", "verbose", "set to: 0 (results only), 1 (confusion matrix), 2 (debug)",0);
+ Optionpk<short> verbose_opt("v", "verbose", "set to: 0 (results only), 1 (confusion matrix), 2 (debug)",0,2);
+
+ tlayer_opt.setHide(1);
+ label_opt.setHide(1);
+ balance_opt.setHide(1);
+ random_opt.setHide(1);
+ minSize_opt.setHide(1);
+ band_opt.setHide(1);
+ bstart_opt.setHide(1);
+ bend_opt.setHide(1);
+ offset_opt.setHide(1);
+ scale_opt.setHide(1);
+ aggreg_opt.setHide(1);
+ // priors_opt.setHide(1);
+ selector_opt.setHide(1);
+ epsilon_cost_opt.setHide(1);
+ cv_opt.setHide(1);
+ classname_opt.setHide(1);
+ classvalue_opt.setHide(1);
+ nneuron_opt.setHide(1);
+ connection_opt.setHide(1);
+ weights_opt.setHide(1);
+ learning_opt.setHide(1);
+ maxit_opt.setHide(1);
bool doProcess;//stop process when program was invoked with help option (-h --help)
try{
@@ -232,6 +255,9 @@ int main(int argc, char *argv[])
exit(0);
}
if(!doProcess){
+ cout << endl;
+ cout << "Usage: pkfsann -t training -n number" << 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
}
diff --git a/src/apps/pkfssvm.cc b/src/apps/pkfssvm.cc
index 83caa76..0ca3ab8 100644
--- a/src/apps/pkfssvm.cc
+++ b/src/apps/pkfssvm.cc
@@ -1,5 +1,5 @@
/**********************************************************************
-pkfssvm.cc: feature selection for svm classifier
+pkfssvm.cc: feature selection for support vector machine classifier pksvm
Copyright (C) 2008-2014 Pieter Kempeneers
This file is part of pktools
@@ -228,7 +228,35 @@ int main(int argc, char *argv[])
Optionpk<unsigned short> cv_opt("cv", "cv", "n-fold cross validation mode",2);
Optionpk<string> classname_opt("c", "class", "list of class names.");
Optionpk<short> classvalue_opt("r", "reclass", "list of class values (use same order as in classname opt.");
- Optionpk<short> verbose_opt("v", "verbose", "set to: 0 (results only), 1 (confusion matrix), 2 (debug)",0);
+ Optionpk<short> verbose_opt("v", "verbose", "set to: 0 (results only), 1 (confusion matrix), 2 (debug)",0,2);
+
+ tlayer_opt.setHide(1);
+ label_opt.setHide(1);
+ balance_opt.setHide(1);
+ random_opt.setHide(1);
+ minSize_opt.setHide(1);
+ band_opt.setHide(1);
+ bstart_opt.setHide(1);
+ bend_opt.setHide(1);
+ offset_opt.setHide(1);
+ scale_opt.setHide(1);
+ svm_type_opt.setHide(1);
+ kernel_type_opt.setHide(1);
+ kernel_degree_opt.setHide(1);
+ gamma_opt.setHide(1);
+ coef0_opt.setHide(1);
+ ccost_opt.setHide(1);
+ nu_opt.setHide(1);
+ epsilon_loss_opt.setHide(1);
+ cache_opt.setHide(1);
+ epsilon_tol_opt.setHide(1);
+ shrinking_opt.setHide(1);
+ prob_est_opt.setHide(1);
+ selector_opt.setHide(1);
+ epsilon_cost_opt.setHide(1);
+ cv_opt.setHide(1);
+ classname_opt.setHide(1);
+ classvalue_opt.setHide(1);
bool doProcess;//stop process when program was invoked with help option (-h --help)
try{
@@ -245,7 +273,6 @@ int main(int argc, char *argv[])
bend_opt.retrieveOption(argc,argv);
offset_opt.retrieveOption(argc,argv);
scale_opt.retrieveOption(argc,argv);
- // priors_opt.retrieveOption(argc,argv);
svm_type_opt.retrieveOption(argc,argv);
kernel_type_opt.retrieveOption(argc,argv);
kernel_degree_opt.retrieveOption(argc,argv);
@@ -270,6 +297,9 @@ int main(int argc, char *argv[])
exit(0);
}
if(!doProcess){
+ cout << endl;
+ cout << "Usage: pkfssvm -t training -n number" << 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
}
diff --git a/src/apps/pkgetmask.cc b/src/apps/pkgetmask.cc
index aa9004e..8ba8b9d 100644
--- a/src/apps/pkgetmask.cc
+++ b/src/apps/pkgetmask.cc
@@ -37,18 +37,25 @@ 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> colorTable_opt("ct", "ct", "color table (file with 5 columns: id R G B ALFA (0: transparent, 255: solid)");
- Optionpk<short> verbose_opt("v", "verbose", "verbose", 0);
+ Optionpk<short> verbose_opt("v", "verbose", "verbose", 0,2);
+
+ band_opt.setHide(1);
+ operator_opt.setHide(1);
+ otype_opt.setHide(1);
+ oformat_opt.setHide(1);
+ option_opt.setHide(1);
+ colorTable_opt.setHide(1);
bool doProcess;//stop process when program was invoked with help option (-h --help)
try{
doProcess=input_opt.retrieveOption(argc,argv);
- band_opt.retrieveOption(argc,argv);
+ output_opt.retrieveOption(argc,argv);
min_opt.retrieveOption(argc,argv);
max_opt.retrieveOption(argc,argv);
- operator_opt.retrieveOption(argc,argv);
data_opt.retrieveOption(argc,argv);
nodata_opt.retrieveOption(argc,argv);
- output_opt.retrieveOption(argc,argv);
+ band_opt.retrieveOption(argc,argv);
+ operator_opt.retrieveOption(argc,argv);
otype_opt.retrieveOption(argc,argv);
oformat_opt.retrieveOption(argc,argv);
option_opt.retrieveOption(argc,argv);
@@ -60,6 +67,9 @@ int main(int argc,char **argv) {
exit(0);
}
if(!doProcess){
+ cout << endl;
+ cout << "Usage: pkgetmask -i input -o output" << 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
}
diff --git a/src/apps/pkinfo.cc b/src/apps/pkinfo.cc
index 2978eed..e244675 100644
--- a/src/apps/pkinfo.cc
+++ b/src/apps/pkinfo.cc
@@ -115,6 +115,9 @@ int main(int argc, char *argv[])
exit(0);
}
if(!doProcess){
+ cout << endl;
+ cout << "Usage: pkinfo -i input [options]" << 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
}
diff --git a/src/apps/pkkalman.cc b/src/apps/pkkalman.cc
index be1de0f..96885db 100644
--- a/src/apps/pkkalman.cc
+++ b/src/apps/pkkalman.cc
@@ -50,10 +50,10 @@ int main(int argc,char **argv) {
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> deltaObs_opt("dobs", "deltaobs", "Threshold for relative difference (in percentage) in observation and model pixel");
+ 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<double> regSensor_opt("rs", "regsensor", "Relative Weight for regression in sensor 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<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);
@@ -87,7 +87,7 @@ int main(int argc,char **argv) {
deltaObs_opt.retrieveOption(argc,argv);
uncertNodata_opt.retrieveOption(argc,argv);
// regTime_opt.retrieveOption(argc,argv);
- // regSensor_opt.retrieveOption(argc,argv);
+ regSensor_opt.retrieveOption(argc,argv);
down_opt.retrieveOption(argc,argv);
threshold_opt.retrieveOption(argc,argv);
minreg_opt.retrieveOption(argc,argv);
@@ -117,6 +117,12 @@ 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;
@@ -205,13 +211,13 @@ int main(int argc,char **argv) {
imgreg.setDown(down_opt[0]);
imgreg.setThreshold(threshold_opt[0]);
- double c0modGlobal=0;//time regression coefficient c0 calculated on entire image
- double c1modGlobal=0;//time regression coefficient c1 calculated on entire image
- double c0mod=0;//time regression coefficient c0 calculated on local window
- double c1mod=0;//time regression coefficient c1 calculated on local window
+ double c0modGlobal=0;//time regression coefficient c0 (offset) calculated on entire image
+ double c1modGlobal=1;//time regression coefficient c1 (scale) calculated on entire image
+ double c0mod=0;//time regression coefficient c0 (offset) calculated on local window
+ double c1mod=1;//time regression coefficient c1 (scale) calculated on local window
- double c0obs=0;
- double c1obs=0;
+ double c0obs=0;//offset
+ double c1obs=1;//scale
double errObs=uncertNodata_opt[0];//start with high initial value in case we do not have first observation at time 0
vector<int> relobsindex;
@@ -233,9 +239,8 @@ int main(int argc,char **argv) {
if(find(direction_opt.begin(),direction_opt.end(),"forward")!=direction_opt.end()){
///////////////////////////// forward model /////////////////////////
+ cout << "Running forward model" << endl;
obsindex=0;
- if(verbose_opt[0])
- cout << "Running forward model" << endl;
//initialization
string output;
if(outputfw_opt.size()==model_opt.size())
@@ -302,8 +307,15 @@ int main(int argc,char **argv) {
for(int icol=0;icol<ncol;++icol){
imgWriterEst.image2geo(icol,irow,x,y);
imgReaderModel1.geo2image(x,y,modCol,modRow);
- estWriteBuffer[icol]=estReadBuffer[modCol];
- uncertWriteBuffer[icol]=uncertModel_opt[0]*stdDev;
+ double modValue=estReadBuffer[modCol];
+ if(imgReaderModel1.isNoData(modValue)){
+ estWriteBuffer[icol]=obsnodata_opt[0];
+ uncertWriteBuffer[icol]=uncertNodata_opt[0];
+ }
+ else{
+ estWriteBuffer[icol]=modValue;
+ uncertWriteBuffer[icol]=uncertModel_opt[0]*stdDev;
+ }
}
imgWriterEst.writeData(estWriteBuffer,GDT_Float64,irow,0);
imgWriterEst.writeData(uncertWriteBuffer,GDT_Float64,irow,1);
@@ -324,51 +336,107 @@ int main(int argc,char **argv) {
imgReaderObs.setNoData(obsnodata_opt);
imgReaderObs.setOffset(obsoffset_opt[0]);
imgReaderObs.setScale(obsscale_opt[0]);
+
+ if(regSensor_opt[0])
+ errObs=imgreg.getRMSE(imgReaderModel1,imgReaderObs,c0obs,c1obs,verbose_opt[0]);
+ else{
+ c0obs=0;
+ c1obs=1;
+ errObs=0;
+ }
+
for(int irow=0;irow<nrow;++irow){
vector<double> estReadBuffer;
imgWriterEst.image2geo(0,irow,x,y);
imgReaderModel1.geo2image(x,y,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> uncertObsBuffer;
- imgReaderObs.readData(estWriteBuffer,GDT_Float64,irow,0);
+ vector<double> uncertObsLineBuffer;
+ // imgReaderObs.readData(estWriteBuffer,GDT_Float64,irow,0);
+ imgReaderObs.readData(obsLineBuffer,GDT_Float64,irow,0);
+
if(imgReaderObs.nrOfBand()>1)
- imgReaderObs.readData(uncertObsBuffer,GDT_Float64,irow,1);
+ imgReaderObs.readData(uncertObsLineBuffer,GDT_Float64,irow,1);
for(int icol=0;icol<ncol;++icol){
- if(imgReaderObs.isNoData(estWriteBuffer[icol])){
- imgWriterEst.image2geo(icol,irow,x,y);
- imgReaderModel1.geo2image(x,y,modCol,modRow);
- assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
- if(imgReaderModel1.isNoData(estReadBuffer[modCol]))//if both obs and model are no-data, set obs to nodata
+ imgWriterEst.image2geo(icol,irow,x,y);
+ imgReaderModel1.geo2image(x,y,modCol,modRow);
+ assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
+ double modValue=estReadBuffer[modCol];
+ if(imgReaderModel1.isNoData(modValue)){//model is nodata: retain observation
+ estWriteBuffer[icol]=obsLineBuffer[icol];
+ if(imgReaderObs.isNoData(obsLineBuffer[icol])){
estWriteBuffer[icol]=obsnodata_opt[0];
+ uncertWriteBuffer[icol]=uncertNodata_opt[0];
+ }
+ else if(uncertObsLineBuffer.size()>icol)
+ uncertWriteBuffer[icol]=uncertObsLineBuffer[icol];
else
- estWriteBuffer[icol]=estReadBuffer[modCol];
- uncertWriteBuffer[icol]=uncertModel_opt[0]*stdDev;
+ uncertWriteBuffer[icol]=uncertObs_opt[0];
}
- else{
+ else{//model is valid: calculate estimate from model
+ double errMod=uncertModel_opt[0]*stdDev;
+ double certNorm=(errMod*errMod+errObs*errObs);
+ double certMod=errObs*errObs/certNorm;
+ double certObs=errMod*errMod/certNorm;
+ double regTime=0;
+ double regSensor=(c0obs+c1obs*modValue)*certObs;
+ estWriteBuffer[icol]=regTime+regSensor;
+ double totalUncertainty=0;
+ if(errMod<eps_opt[0])
+ totalUncertainty=errObs;
+ else if(errObs<eps_opt[0])
+ totalUncertainty=errMod;
+ else{
+ totalUncertainty=1.0/errMod/errMod+1/errObs/errObs;
+ totalUncertainty=sqrt(1.0/totalUncertainty);
+ }
+ uncertWriteBuffer[icol]=totalUncertainty;//in case observation is not valid
+ }
+ //todo: check with Fernando? (here uncertainty only relates to modeled estimate. In case of observation update, we include uncertainty of observation
+ // uncertWriteBuffer[icol]=totalUncertainty+uncertReadBuffer[icol];
+ //observation update if observation is valid
+
+ if(!imgReaderObs.isNoData(obsLineBuffer[icol])){
+ double kalmanGain=1;
double uncertObs=uncertObs_opt[0];
- if(uncertObsBuffer.size()>icol)
- uncertObs=uncertObsBuffer[icol];
- if(uncertModel_opt[0]*stdDev+uncertObs>eps_opt[0]){
- imgWriterEst.image2geo(icol,irow,x,y);
- imgReaderModel1.geo2image(x,y,modCol,modRow);
- assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
- // double noemer=uncertObs*uncertObs+stdDev*stdDev;//todo: multiply stdDev with uncertModel_opt[0]
- // estWriteBuffer[icol]*=uncertModel_opt[0]*stdDev*stdDev/noemer;
- // estWriteBuffer[icol]+=uncertModel_opt[0]*uncertObs*uncertObs/noemer;//todo:check! error?
- double kalmanGain=1;
- if(!imgReaderModel1.isNoData(estReadBuffer[modCol])){
- //if model is no-data, retain observation value
- kalmanGain=uncertModel_opt[0]*stdDev/(uncertModel_opt[0]*stdDev+uncertObs);
- estWriteBuffer[icol]=estReadBuffer[modCol]+kalmanGain*(estWriteBuffer[icol]-estReadBuffer[modCol]);
+ bool doUpdate=true;
+ if(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;
+ int minRow=(irow>down_opt[0]/2) ? irow-down_opt[0]/2 : 0;
+ int maxRow=(irow+down_opt[0]/2<imgReaderObs.nrOfRow()) ? irow+down_opt[0]/2 : imgReaderObs.nrOfRow()-1;
+
+ imgReaderObs.readDataBlock(obsWindowBuffer,GDT_Float64,minCol,maxCol,minRow,maxRow,0);
+ statfactory::StatFactory statobs;
+ 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;
}
- uncertWriteBuffer[icol]=uncertModel_opt[0]*stdDev*(1-kalmanGain);
}
- else{
- //no need to fill write buffer (already done in imgReaderObs.readData
- uncertWriteBuffer[icol]=uncertObs;
+ if(doUpdate){
+ if(uncertObsLineBuffer.size()>icol)
+ uncertObs=uncertObsLineBuffer[icol];
+ 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);
}
}
}
@@ -442,7 +510,13 @@ int main(int argc,char **argv) {
//calculate regression between model and observation
if(verbose_opt[0])
cout << "Calculating regression for " << imgReaderModel2.getFileName() << " " << imgReaderObs.getFileName() << endl;
- errObs=imgreg.getRMSE(imgReaderModel2,imgReaderObs,c0obs,c1obs,verbose_opt[0]);
+ if(regSensor_opt[0])
+ errObs=imgreg.getRMSE(imgReaderModel2,imgReaderObs,c0obs,c1obs,verbose_opt[0]);
+ else{
+ c0obs=0;
+ c1obs=1;
+ errObs=0;
+ }
if(verbose_opt[0])
cout << "c0obs, c1obs: " << c0obs << ", " << c1obs << endl;
}
@@ -460,6 +534,7 @@ int main(int argc,char **argv) {
imgReaderEst.setOffset(obsoffset_opt[0]);
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
vector<double> model1LineBuffer;
@@ -473,9 +548,17 @@ int main(int argc,char **argv) {
vector<double> estWriteBuffer(ncol);
vector<double> uncertWriteBuffer(ncol);
+ //initialize obsLineVector
+ assert(down_opt[0]%2);//window size must be odd
+ for(int iline=-down_opt[0]/2+1;iline<down_opt[0]/2+1;++iline){
+ if(iline<0)//replicate line 0
+ imgReaderObs.readData(obsLineVector[iline+down_opt[0]/2],GDT_Float64,0,0);
+ else
+ imgReaderObs.readData(obsLineVector[iline+down_opt[0]/2],GDT_Float64,iline,0);
+ }
for(int irow=0;irow<imgWriterEst.nrOfRow();++irow){
assert(irow<imgReaderEst.nrOfRow());
- //not needed here, because we read entire window for each pixel...
+ //do not read from imgReaderObs, because we read entire window for each pixel...
imgReaderEst.readData(estReadBuffer,GDT_Float64,irow,0);
imgReaderEst.readData(uncertReadBuffer,GDT_Float64,irow,1);
//read model2 in case current estimate is nodata
@@ -489,7 +572,12 @@ int main(int argc,char **argv) {
imgReaderModel1.readData(model1LineBuffer,GDT_Float64,modRow,0);
if(update){
- imgReaderObs.readData(obsLineBuffer,GDT_Float64,irow,0);
+ int maxRow=(irow+down_opt[0]/2<imgReaderEst.nrOfRow()) ? irow+down_opt[0]/2 : imgReaderEst.nrOfRow()-1;
+ obsLineVector.erase(obsLineVector.begin());
+ imgReaderObs.readData(obsLineBuffer,GDT_Float64,maxRow,0);
+ obsLineVector.push_back(obsLineBuffer);
+ obsLineBuffer=obsLineVector[down_opt[0]/2];
+ // imgReaderObs.readData(obsLineBuffer,GDT_Float64,irow,0);
if(imgReaderObs.nrOfBand()>1)
imgReaderObs.readData(uncertObsLineBuffer,GDT_Float64,irow,1);
}
@@ -498,11 +586,17 @@ int main(int argc,char **argv) {
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;
int maxRow=(irow+down_opt[0]/2<imgReaderEst.nrOfRow()) ? irow+down_opt[0]/2 : imgReaderEst.nrOfRow()-1;
- if(update)
- imgReaderObs.readDataBlock(obsWindowBuffer,GDT_Float64,minCol,maxCol,minRow,maxRow,0);
- // imgReaderEst.readDataBlock(estWindowBuffer,GDT_Float64,minCol,maxCol,minRow,maxRow,0);
+ if(update){
+ obsWindowBuffer.clear();
+ for(int iline=0;iline<obsLineVector.size();++iline){
+ for(int isample=minCol;isample<=maxCol;++isample){
+ assert(isample<obsLineVector[iline].size());
+ obsWindowBuffer.push_back(obsLineVector[iline][isample]);
+ }
+ }
+ // imgReaderObs.readDataBlock(obsWindowBuffer,GDT_Float64,minCol,maxCol,minRow,maxRow,0);
+ }
double estValue=estReadBuffer[icol];
- // double estValue=estWindowBuffer[estWindowBuffer.size()/2];
double estMeanValue=0;//stat.mean(estWindowBuffer);
double nvalid=0;
//time update
@@ -510,22 +604,7 @@ int main(int argc,char **argv) {
imgReaderModel2.geo2image(x,y,modCol,modRow);
assert(modRow>=0&&modRow<imgReaderModel2.nrOfRow());
double modValue=model2LineBuffer[modCol];
- bool estNodata=imgReaderEst.isNoData(estValue);
- //todo: debug checkDiff_opt
- // if(checkDiff_opt[0]){
- // vector<double>::iterator itwin=estWindowBuffer.begin();
- // while(itwin!=estWindowBuffer.end()){
- // if(!imgReaderEst.isNoData(*itwin)){
- // estMeanValue+=*itwin;
- // ++nvalid;
- // }
- // ++itwin;
- // }
- // estMeanValue/=nvalid;
- // double difference=(estMeanValue-c0obs)/c1obs-model1LineBuffer[modCol];
- // difference*=difference;
- // estNodata=estNodata||(sqrt(difference)>uncertModel_opt[0]*stdDev);
- // }
+ bool estNodata=imgReaderEst.isNoData(estValue);//validity of current estimate
if(estNodata){
//we have not found any valid data yet, better here to take the current model value if valid
if(imgReaderModel2.isNoData(modValue)){//if both estimate and model are no-data, set obs to nodata
@@ -567,12 +646,6 @@ int main(int argc,char **argv) {
bool modNodata=false;
modNodata=modNodata||imgReaderModel1.isNoData(*it1);
modNodata=modNodata||imgReaderModel2.isNoData(*it2);
- //todo: debug checkDiff_opt
- // if(checkDiff_opt[0]){
- // double difference=(estMeanValue-c0obs)/c1obs-*it1;
- // difference*=difference;
- // modNodata=modNodata||(sqrt(difference)>uncertModel_opt[0]*stdDev);
- // }
if(modNodata){
model1buffer.erase(it1);
model2buffer.erase(it2);
@@ -610,35 +683,40 @@ int main(int argc,char **argv) {
totalUncertainty=1.0/errMod/errMod+1/errObs/errObs;
totalUncertainty=sqrt(1.0/totalUncertainty);
}
- uncertWriteBuffer[icol]=totalUncertainty+uncertReadBuffer[icol];
+ // uncertWriteBuffer[icol]=totalUncertainty+uncertReadBuffer[icol];
+ uncertWriteBuffer[icol]=totalUncertainty;
}
//observation update
if(update&&!imgReaderObs.isNoData(obsLineBuffer[icol])){
+ double kalmanGain=1;
+ double uncertObs=uncertObs_opt[0];
bool doUpdate=true;
if(deltaObs_opt.size()){
statfactory::StatFactory statobs;
statobs.setNoDataValues(obsnodata_opt);
double obsMeanValue=statobs.mean(obsWindowBuffer);
double difference=(obsMeanValue-c0obs)/c1obs-modValue;
- difference/=modValue;//make relative difference
- if(deltaObs_opt[0]<0)
- doUpdate=((100*difference)<deltaObs_opt[0]);
- else{
- difference*=difference;
- doUpdate=(100*sqrt(difference)<deltaObs_opt[0]);
+ 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(doUpdate){
- double kalmanGain=1;
- double uncertObs=uncertObs_opt[0];
if(uncertObsLineBuffer.size()>icol)
uncertObs=uncertObsLineBuffer[icol];
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);
}
+ else
+ kalmanGain=0;
+ estWriteBuffer[icol]+=kalmanGain*(obsLineBuffer[icol]-estWriteBuffer[icol]);
+ uncertWriteBuffer[icol]*=(1-kalmanGain);
}
}
imgWriterEst.writeData(estWriteBuffer,GDT_Float64,irow,0);
@@ -660,9 +738,8 @@ int main(int argc,char **argv) {
}
if(find(direction_opt.begin(),direction_opt.end(),"backward")!=direction_opt.end()){
///////////////////////////// backward model /////////////////////////
+ cout << "Running backward model" << endl;
obsindex=relobsindex.size()-1;
- if(verbose_opt[0])
- cout << "Running backward model" << endl;
//initialization
string output;
if(outputbw_opt.size()==model_opt.size())
@@ -751,51 +828,107 @@ int main(int argc,char **argv) {
imgReaderObs.setNoData(obsnodata_opt);
imgReaderObs.setOffset(obsoffset_opt[0]);
imgReaderObs.setScale(obsscale_opt[0]);
+
+ if(regSensor_opt[0])
+ errObs=imgreg.getRMSE(imgReaderModel1,imgReaderObs,c0obs,c1obs,verbose_opt[0]);
+ else{
+ c0obs=0;
+ c1obs=1;
+ errObs=0;
+ }
+
for(int irow=0;irow<nrow;++irow){
vector<double> estReadBuffer;
imgWriterEst.image2geo(0,irow,x,y);
imgReaderModel1.geo2image(x,y,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;
- imgReaderObs.readData(estWriteBuffer,GDT_Float64,irow,0);
+ // 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);
for(int icol=0;icol<imgWriterEst.nrOfCol();++icol){
- if(imgReaderObs.isNoData(estWriteBuffer[icol])){
- imgWriterEst.image2geo(icol,irow,x,y);
- imgReaderModel1.geo2image(x,y,modCol,modRow);
- assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
- if(imgReaderModel1.isNoData(estReadBuffer[modCol]))//if both obs and model are no-data, set obs to nodata
+ imgWriterEst.image2geo(icol,irow,x,y);
+ imgReaderModel1.geo2image(x,y,modCol,modRow);
+ assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
+ double modValue=estReadBuffer[modCol];
+ if(imgReaderModel1.isNoData(modValue)){//model is nodata: retain observation
+ estWriteBuffer[icol]=obsLineBuffer[icol];
+ if(imgReaderObs.isNoData(obsLineBuffer[icol])){
estWriteBuffer[icol]=obsnodata_opt[0];
+ uncertWriteBuffer[icol]=uncertNodata_opt[0];
+ }
+ else if(uncertObsLineBuffer.size()>icol)
+ uncertWriteBuffer[icol]=uncertObsLineBuffer[icol];
else
- estWriteBuffer[icol]=estReadBuffer[modCol];
- uncertWriteBuffer[icol]=uncertModel_opt[0]*stdDev;
+ uncertWriteBuffer[icol]=uncertObs_opt[0];
}
- else{
+ else{//model is valid: calculate estimate from model
+ double errMod=uncertModel_opt[0]*stdDev;
+ double certNorm=(errMod*errMod+errObs*errObs);
+ double certMod=errObs*errObs/certNorm;
+ double certObs=errMod*errMod/certNorm;
+ double regTime=0;
+ double regSensor=(c0obs+c1obs*modValue)*certObs;
+ estWriteBuffer[icol]=regTime+regSensor;
+ double totalUncertainty=0;
+ if(errMod<eps_opt[0])
+ totalUncertainty=errObs;
+ else if(errObs<eps_opt[0])
+ totalUncertainty=errMod;
+ else{
+ totalUncertainty=1.0/errMod/errMod+1/errObs/errObs;
+ totalUncertainty=sqrt(1.0/totalUncertainty);
+ }
+ uncertWriteBuffer[icol]=totalUncertainty;//in case observation is not valid
+ }
+ //todo: check with Fernando? (here uncertainty only relates to modeled estimate. In case of observation update, we include uncertainty of observation
+ // uncertWriteBuffer[icol]=totalUncertainty+uncertReadBuffer[icol];
+ //observation update if observation is valid
+
+ if(!imgReaderObs.isNoData(obsLineBuffer[icol])){
+ double kalmanGain=1;
+ bool doUpdate=true;
double uncertObs=uncertObs_opt[0];
- if(uncertObsLineBuffer.size()>icol)
- uncertObs=uncertObsLineBuffer[icol];
- if(uncertModel_opt[0]*stdDev+uncertObs>eps_opt[0]){
- imgWriterEst.image2geo(icol,irow,x,y);
- imgReaderModel1.geo2image(x,y,modCol,modRow);
- assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
- // double noemer=uncertObs*uncertObs+stdDev*stdDev;//todo: multiply stdDev with uncertModel_opt[0]
- // estWriteBuffer[icol]*=uncertModel_opt[0]*stdDev*stdDev/noemer;
- // estWriteBuffer[icol]+=uncertModel_opt[0]*uncertObs*uncertObs/noemer;//todo:check! error?
- double kalmanGain=1;
- if(!imgReaderModel1.isNoData(estReadBuffer[modCol])){
- //if model is no-data, retain observation value
- kalmanGain=uncertModel_opt[0]*stdDev/(uncertModel_opt[0]*stdDev+uncertObs);
- estWriteBuffer[icol]=estReadBuffer[modCol]+kalmanGain*(estWriteBuffer[icol]-estReadBuffer[modCol]);
+ if(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;
+ int minRow=(irow>down_opt[0]/2) ? irow-down_opt[0]/2 : 0;
+ int maxRow=(irow+down_opt[0]/2<imgReaderObs.nrOfRow()) ? irow+down_opt[0]/2 : imgReaderObs.nrOfRow()-1;
+
+ imgReaderObs.readDataBlock(obsWindowBuffer,GDT_Float64,minCol,maxCol,minRow,maxRow,0);
+ statfactory::StatFactory statobs;
+ 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;
}
- uncertWriteBuffer[icol]=uncertModel_opt[0]*stdDev*(1-kalmanGain);
}
- else{
- //no need to fill write buffer (already done in imgReaderObs.readData
- uncertWriteBuffer[icol]=uncertObs;
+ if(doUpdate){
+ if(uncertObsLineBuffer.size()>icol)
+ uncertObs=uncertObsLineBuffer[icol];
+ 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);
}
}
}
@@ -869,7 +1002,13 @@ int main(int argc,char **argv) {
//calculate regression between model and observation
if(verbose_opt[0])
cout << "Calculating regression for " << imgReaderModel2.getFileName() << " " << imgReaderObs.getFileName() << endl;
- errObs=imgreg.getRMSE(imgReaderModel2,imgReaderObs,c0obs,c1obs,verbose_opt[0]);
+ if(regSensor_opt[0])
+ errObs=imgreg.getRMSE(imgReaderModel2,imgReaderObs,c0obs,c1obs,verbose_opt[0]);
+ else{
+ c0obs=0;
+ c1obs=1;
+ errObs=0;
+ }
if(verbose_opt[0])
cout << "c0obs, c1obs: " << c0obs << ", " << c1obs << endl;
}
@@ -887,6 +1026,7 @@ int main(int argc,char **argv) {
imgReaderEst.setOffset(obsoffset_opt[0]);
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
vector<double> model1LineBuffer;
@@ -900,9 +1040,17 @@ int main(int argc,char **argv) {
vector<double> estWriteBuffer(ncol);
vector<double> uncertWriteBuffer(ncol);
+ //initialize obsLineVector
+ assert(down_opt[0]%2);//window size must be odd
+ for(int iline=-down_opt[0]/2+1;iline<down_opt[0]/2+1;++iline){
+ if(iline<0)//replicate line 0
+ imgReaderObs.readData(obsLineVector[iline+down_opt[0]/2],GDT_Float64,0,0);
+ else
+ imgReaderObs.readData(obsLineVector[iline+down_opt[0]/2],GDT_Float64,iline,0);
+ }
for(int irow=0;irow<imgWriterEst.nrOfRow();++irow){
assert(irow<imgReaderEst.nrOfRow());
- //not needed here, because we read entire window for each pixel...
+ //do not read from imgReaderObs, because we read entire window for each pixel...
imgReaderEst.readData(estReadBuffer,GDT_Float64,irow,0);
imgReaderEst.readData(uncertReadBuffer,GDT_Float64,irow,1);
//read model2 in case current estimate is nodata
@@ -916,7 +1064,12 @@ int main(int argc,char **argv) {
imgReaderModel1.readData(model1LineBuffer,GDT_Float64,modRow,0);
if(update){
- imgReaderObs.readData(obsLineBuffer,GDT_Float64,irow,0);
+ int maxRow=(irow+down_opt[0]/2<imgReaderEst.nrOfRow()) ? irow+down_opt[0]/2 : imgReaderEst.nrOfRow()-1;
+ obsLineVector.erase(obsLineVector.begin());
+ imgReaderObs.readData(obsLineBuffer,GDT_Float64,maxRow,0);
+ obsLineVector.push_back(obsLineBuffer);
+ obsLineBuffer=obsLineVector[down_opt[0]/2];
+ // imgReaderObs.readData(obsLineBuffer,GDT_Float64,irow,0);
if(imgReaderObs.nrOfBand()>1)
imgReaderObs.readData(uncertObsLineBuffer,GDT_Float64,irow,1);
}
@@ -925,11 +1078,17 @@ int main(int argc,char **argv) {
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;
int maxRow=(irow+down_opt[0]/2<imgReaderEst.nrOfRow()) ? irow+down_opt[0]/2 : imgReaderEst.nrOfRow()-1;
- if(update)
- imgReaderObs.readDataBlock(obsWindowBuffer,GDT_Float64,minCol,maxCol,minRow,maxRow,0);
- // imgReaderEst.readDataBlock(estWindowBuffer,GDT_Float64,minCol,maxCol,minRow,maxRow,0);
+ if(update){
+ obsWindowBuffer.clear();
+ for(int iline=0;iline<obsLineVector.size();++iline){
+ for(int isample=minCol;isample<=maxCol;++isample){
+ assert(isample<obsLineVector[iline].size());
+ obsWindowBuffer.push_back(obsLineVector[iline][isample]);
+ }
+ }
+ // imgReaderObs.readDataBlock(obsWindowBuffer,GDT_Float64,minCol,maxCol,minRow,maxRow,0);
+ }
double estValue=estReadBuffer[icol];
- // double estValue=estWindowBuffer[estWindowBuffer.size()/2];
double estMeanValue=0;//stat.mean(estWindowBuffer);
double nvalid=0;
//time update
@@ -938,21 +1097,6 @@ int main(int argc,char **argv) {
assert(modRow>=0&&modRow<imgReaderModel2.nrOfRow());
double modValue=model2LineBuffer[modCol];
bool estNodata=imgReaderEst.isNoData(estValue);
- //todo: debug checkDiff_opt
- // if(checkDiff_opt[0]){
- // vector<double>::iterator itwin=estWindowBuffer.begin();
- // while(itwin!=estWindowBuffer.end()){
- // if(!imgReaderEst.isNoData(*itwin)){
- // estMeanValue+=*itwin;
- // ++nvalid;
- // }
- // ++itwin;
- // }
- // estMeanValue/=nvalid;
- // double difference=(estMeanValue-c0obs)/c1obs-model1LineBuffer[modCol];
- // difference*=difference;
- // estNodata=estNodata||(sqrt(difference)>uncertModel_opt[0]*stdDev);
- // }
if(estNodata){
//we have not found any valid data yet, better here to take the current model value if valid
if(imgReaderModel2.isNoData(modValue)){//if both estimate and model are no-data, set obs to nodata
@@ -994,12 +1138,6 @@ int main(int argc,char **argv) {
bool modNodata=false;
modNodata=modNodata||imgReaderModel1.isNoData(*it1);
modNodata=modNodata||imgReaderModel2.isNoData(*it2);
- //todo: debug checkDiff_opt
- // if(checkDiff_opt[0]){
- // double difference=(estMeanValue-c0obs)/c1obs-*it1;
- // difference*=difference;
- // modNodata=modNodata||(sqrt(difference)>uncertModel_opt[0]*stdDev);
- // }
if(modNodata){
model1buffer.erase(it1);
model2buffer.erase(it2);
@@ -1009,7 +1147,7 @@ int main(int argc,char **argv) {
++it2;
}
}
- if(model1buffer.size()>minreg_opt[0]&&model2buffer.size()>minreg_opt[5])
+ if(model1buffer.size()>minreg_opt[0]&&model2buffer.size()>minreg_opt[0])
errMod=stat.linear_regression_err(model1buffer,model2buffer,c0mod,c1mod);
else{//use global regression...
c0mod=c0modGlobal;
@@ -1037,37 +1175,40 @@ int main(int argc,char **argv) {
totalUncertainty=1.0/errMod/errMod+1/errObs/errObs;
totalUncertainty=sqrt(1.0/totalUncertainty);
}
- uncertWriteBuffer[icol]=totalUncertainty+uncertReadBuffer[icol];
+ // uncertWriteBuffer[icol]=totalUncertainty+uncertReadBuffer[icol];
+ uncertWriteBuffer[icol]=totalUncertainty;
}
//observation update
if(update&&!imgReaderObs.isNoData(obsLineBuffer[icol])){
+ double kalmanGain=1;
+ double uncertObs=uncertObs_opt[0];
bool doUpdate=true;
if(deltaObs_opt.size()){
statfactory::StatFactory statobs;
statobs.setNoDataValues(obsnodata_opt);
double obsMeanValue=statobs.mean(obsWindowBuffer);
double difference=(obsMeanValue-c0obs)/c1obs-modValue;
- difference*=difference;
- difference/=modValue;//make relative difference
- difference*=difference;
- if(deltaObs_opt[0]<0)
- doUpdate=((100*difference)<deltaObs_opt[0]);
- else{
- difference*=difference;
- doUpdate=(100*sqrt(difference)<deltaObs_opt[0]);
+ 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(doUpdate){
- double kalmanGain=1;
- double uncertObs=uncertObs_opt[0];
if(uncertObsLineBuffer.size()>icol)
uncertObs=uncertObsLineBuffer[icol];
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);
}
+ else
+ kalmanGain=0;
+ estWriteBuffer[icol]+=kalmanGain*(obsLineBuffer[icol]-estWriteBuffer[icol]);
+ uncertWriteBuffer[icol]*=(1-kalmanGain);
}
}
imgWriterEst.writeData(estWriteBuffer,GDT_Float64,irow,0);
@@ -1089,9 +1230,8 @@ int main(int argc,char **argv) {
}
if(find(direction_opt.begin(),direction_opt.end(),"smooth")!=direction_opt.end()){
///////////////////////////// smooth model /////////////////////////
+ cout << "Running smooth model" << endl;
obsindex=0;
- if(verbose_opt[0])
- cout << "Running smooth model" << endl;
for(int modindex=0;modindex<model_opt.size();++modindex){
if(verbose_opt[0]){
cout << "processing time " << tmodel_opt[modindex] << endl;
diff --git a/src/apps/pklas2img.cc b/src/apps/pklas2img.cc
index ff6f542..134a385 100644
--- a/src/apps/pklas2img.cc
+++ b/src/apps/pklas2img.cc
@@ -30,13 +30,11 @@ using namespace std;
int main(int argc,char **argv) {
Optionpk<string> input_opt("i", "input", "Input las file");
- Optionpk<short> nodata_opt("nodata", "nodata", "nodata value to put in image if not valid", 0);
Optionpk<string> attribute_opt("n", "name", "names of the attribute to select: intensity, return, nreturn, z", "z");
// Optionpk<bool> disc_opt("circ", "circular", "circular disc kernel for dilation and erosion", false);
// Optionpk<double> maxSlope_opt("s", "maxSlope", "Maximum slope used for morphological filtering", 0.0);
// Optionpk<double> hThreshold_opt("ht", "maxHeight", "initial and maximum height threshold for progressive morphological filtering (e.g., -ht 0.2 -ht 2.5)", 0.2);
// Optionpk<short> maxIter_opt("maxit", "maxit", "Maximum number of iterations in post filter", 5);
- Optionpk<short> nbin_opt("nbin", "nbin", "Number of percentile bins for calculating percentile height value profile (=number of output bands)", 10.0);
Optionpk<unsigned short> returns_opt("ret", "ret", "number(s) of returns to include");
Optionpk<unsigned short> classes_opt("class", "class", "classes to keep: 0 (created, never classified), 1 (unclassified), 2 (ground), 3 (low vegetation), 4 (medium vegetation), 5 (high vegetation), 6 (building), 7 (low point, noise), 8 (model key-point), 9 (water), 10 (reserved), 11 (reserved), 12 (overlap)");
Optionpk<string> composite_opt("comp", "comp", "composite for multiple points in cell (min, max, median, mean, sum, first, last, profile (percentile height values), number (point density)). Last: overwrite cells with latest point", "last");
@@ -52,31 +50,27 @@ int main(int argc,char **argv) {
Optionpk<double> lry_opt("lry", "lry", "Lower right y value bounding box (in geocoordinates if georef is true). 0 is read from input file", 0.0);
Optionpk<string> otype_opt("ot", "otype", "Data type for output image ({Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/CInt16/CInt32/CFloat32/CFloat64}). Empty string: inherit type from input image", "Byte");
Optionpk<string> oformat_opt("of", "oformat", "Output image format (see also gdal_translate). Empty string: inherit from input image", "GTiff");
- Optionpk<string> option_opt("co", "co", "Creation option for output file. Multiple options can be specified.");
Optionpk<double> dx_opt("dx", "dx", "Output resolution in x (in meter)", 1.0);
Optionpk<double> dy_opt("dy", "dy", "Output resolution in y (in meter)", 1.0);
+ Optionpk<short> nbin_opt("nbin", "nbin", "Number of percentile bins for calculating percentile height value profile (=number of output bands)", 10.0);
+ Optionpk<short> nodata_opt("nodata", "nodata", "nodata value to put in image if not valid", 0);
+ Optionpk<string> option_opt("co", "co", "Creation option for output file. Multiple options can be specified.");
Optionpk<string> colorTable_opt("ct", "ct", "color table (file with 5 columns: id R G B ALFA (0: transparent, 255: solid)");
- Optionpk<short> verbose_opt("v", "verbose", "verbose mode", 0);
+ Optionpk<short> verbose_opt("v", "verbose", "verbose mode", 0,2);
+
+ nbin_opt.setHide(1);
+ nodata_opt.setHide(1);
+ option_opt.setHide(1);
+ colorTable_opt.setHide(1);
bool doProcess;//stop process when program was invoked with help option (-h --help)
try{
doProcess=input_opt.retrieveOption(argc,argv);
- // mask_opt.retrieveOption(argc,argv);
- // invalid_opt.retrieveOption(argc,argv);
- nodata_opt.retrieveOption(argc,argv);
attribute_opt.retrieveOption(argc,argv);
- // disc_opt.retrieveOption(argc,argv);
- // maxSlope_opt.retrieveOption(argc,argv);
- // hThreshold_opt.retrieveOption(argc,argv);
- // maxIter_opt.retrieveOption(argc,argv);
- nbin_opt.retrieveOption(argc,argv);
returns_opt.retrieveOption(argc,argv);
classes_opt.retrieveOption(argc,argv);
composite_opt.retrieveOption(argc,argv);
filter_opt.retrieveOption(argc,argv);
- // postFilter_opt.retrieveOption(argc,argv);
- // dimx_opt.retrieveOption(argc,argv);
- // dimy_opt.retrieveOption(argc,argv);
output_opt.retrieveOption(argc,argv);
projection_opt.retrieveOption(argc,argv);
ulx_opt.retrieveOption(argc,argv);
@@ -85,9 +79,11 @@ int main(int argc,char **argv) {
lry_opt.retrieveOption(argc,argv);
otype_opt.retrieveOption(argc,argv);
oformat_opt.retrieveOption(argc,argv);
- option_opt.retrieveOption(argc,argv);
dx_opt.retrieveOption(argc,argv);
dy_opt.retrieveOption(argc,argv);
+ nbin_opt.retrieveOption(argc,argv);
+ nodata_opt.retrieveOption(argc,argv);
+ option_opt.retrieveOption(argc,argv);
colorTable_opt.retrieveOption(argc,argv);
verbose_opt.retrieveOption(argc,argv);
}
@@ -95,7 +91,11 @@ int main(int argc,char **argv) {
std::cout << predefinedString << std::endl;
exit(0);
}
+
if(!doProcess){
+ cout << endl;
+ cout << "pklas2img -i lasfile -o output" << 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
}
diff --git a/src/apps/pkoptsvm.cc b/src/apps/pkoptsvm.cc
index 2463237..a01c249 100644
--- a/src/apps/pkoptsvm.cc
+++ b/src/apps/pkoptsvm.cc
@@ -1,5 +1,5 @@
/**********************************************************************
-pkoptsvm.cc: program to optimize parameters for SVM classification
+pkoptsvm.cc: program to optimize parameters for support vector machine classifier pksvm
Copyright (C) 2008-2014 Pieter Kempeneers
This file is part of pktools
@@ -62,7 +62,7 @@ Optionpk<bool> costfunction_opt("cf", "cf", "use Overall Accuracy instead of kap
Optionpk<unsigned short> cv_opt("cv", "cv", "n-fold cross validation mode",2);
Optionpk<string> classname_opt("c", "class", "list of class names.");
Optionpk<short> classvalue_opt("r", "reclass", "list of class values (use same order as in class opt).");
-Optionpk<short> verbose_opt("v", "verbose", "use 1 to output intermediate results for plotting",0);
+Optionpk<short> verbose_opt("v", "verbose", "use 1 to output intermediate results for plotting",0,2);
double objFunction(const std::vector<double> &x, std::vector<double> &grad, void *my_func_data){
@@ -259,6 +259,10 @@ int main(int argc, char *argv[])
map<short,int> reclassMap;
vector<int> vreclass;
Optionpk<string> training_opt("t", "training", "training vector file. A single vector file contains all training features (must be set as: b0, b1, b2,...) for all classes (class numbers identified by label option).");
+ Optionpk<float> ccost_opt("cc", "ccost", "min and max boundaries the parameter C of C-SVC, epsilon-SVR, and nu-SVR (optional: initial value)",1);
+ Optionpk<float> gamma_opt("g", "gamma", "min max boundaries for gamma in kernel function (optional: initial value)",0);
+ Optionpk<double> stepcc_opt("stepcc","stepcc","multiplicative step for ccost in GRID search",2);
+ Optionpk<double> stepg_opt("stepg","stepg","multiplicative step for gamma in GRID search",2);
Optionpk<string> input_opt("i", "input", "input test vector file");
Optionpk<string> tlayer_opt("tln", "tln", "training layer name(s)");
Optionpk<string> label_opt("label", "label", "identifier for class label in training vector file.","label");
@@ -271,20 +275,49 @@ int main(int argc, char *argv[])
Optionpk<double> bend_opt("e", "end", "end band sequence number (set to 0 to include all bands)", 0);
Optionpk<double> offset_opt("\0", "offset", "offset value for each spectral band input features: refl[band]=(DN[band]-offset[band])/scale[band]", 0.0);
Optionpk<double> scale_opt("\0", "scale", "scale value for each spectral band input features: refl=(DN[band]-offset[band])/scale[band] (use 0 if scale min and max in each band to -1.0 and 1.0)", 0.0);
- Optionpk<float> gamma_opt("g", "gamma", "min max boundaries for gamma in kernel function (optional: initial value)",0);
- Optionpk<float> ccost_opt("cc", "ccost", "min and max boundaries the parameter C of C-SVC, epsilon-SVR, and nu-SVR (optional: initial value)",1);
Optionpk<unsigned int> maxit_opt("maxit","maxit","maximum number of iterations",500);
Optionpk<string> algorithm_opt("a", "algorithm", "GRID, or any optimization algorithm from http://ab-initio.mit.edu/wiki/index.php/NLopt_Algorithms","GRID");
Optionpk<double> tolerance_opt("tol","tolerance","relative tolerance for stopping criterion",0.0001);
- Optionpk<double> step_opt("step","step","multiplicative step for ccost and gamma in GRID search",2);
+ input_opt.setHide(1);
+ tlayer_opt.setHide(1);
+ label_opt.setHide(1);
+ balance_opt.setHide(1);
+ random_opt.setHide(1);
+ minSize_opt.setHide(1);
+ band_opt.setHide(1);
+ bstart_opt.setHide(1);
+ bend_opt.setHide(1);
+ offset_opt.setHide(1);
+ scale_opt.setHide(1);
+ svm_type_opt.setHide(1);
+ kernel_type_opt.setHide(1);
+ kernel_degree_opt.setHide(1);
+ coef0_opt.setHide(1);
+ nu_opt.setHide(1);
+ epsilon_loss_opt.setHide(1);
+ cache_opt.setHide(1);
+ epsilon_tol_opt.setHide(1);
+ shrinking_opt.setHide(1);
+ prob_est_opt.setHide(1);
+ cv_opt.setHide(1);
+ costfunction_opt.setHide(1);
+ maxit_opt.setHide(1);
+ tolerance_opt.setHide(1);
+ algorithm_opt.setHide(1);
+ classname_opt.setHide(1);
+ classvalue_opt.setHide(1);
+
bool doProcess;//stop process when program was invoked with help option (-h --help)
try{
doProcess=training_opt.retrieveOption(argc,argv);
+ ccost_opt.retrieveOption(argc,argv);
+ gamma_opt.retrieveOption(argc,argv);
+ stepcc_opt.retrieveOption(argc,argv);
+ stepg_opt.retrieveOption(argc,argv);
input_opt.retrieveOption(argc,argv);
tlayer_opt.retrieveOption(argc,argv);
label_opt.retrieveOption(argc,argv);
- // reclass_opt.retrieveOption(argc,argv);
balance_opt.retrieveOption(argc,argv);
random_opt.retrieveOption(argc,argv);
minSize_opt.retrieveOption(argc,argv);
@@ -296,9 +329,7 @@ int main(int argc, char *argv[])
svm_type_opt.retrieveOption(argc,argv);
kernel_type_opt.retrieveOption(argc,argv);
kernel_degree_opt.retrieveOption(argc,argv);
- gamma_opt.retrieveOption(argc,argv);
coef0_opt.retrieveOption(argc,argv);
- ccost_opt.retrieveOption(argc,argv);
nu_opt.retrieveOption(argc,argv);
epsilon_loss_opt.retrieveOption(argc,argv);
cache_opt.retrieveOption(argc,argv);
@@ -309,7 +340,6 @@ int main(int argc, char *argv[])
costfunction_opt.retrieveOption(argc,argv);
maxit_opt.retrieveOption(argc,argv);
tolerance_opt.retrieveOption(argc,argv);
- step_opt.retrieveOption(argc,argv);
algorithm_opt.retrieveOption(argc,argv);
classname_opt.retrieveOption(argc,argv);
classvalue_opt.retrieveOption(argc,argv);
@@ -320,6 +350,9 @@ int main(int argc, char *argv[])
exit(0);
}
if(!doProcess){
+ cout << endl;
+ cout << "Usage: pkoptsvm -t training" << 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
}
@@ -625,8 +658,6 @@ int main(int argc, char *argv[])
std::vector<double> x(2);
if(algorithm_opt[0]=="GRID"){
- if(step_opt.size()<2)//[0] for cost, [1] for gamma
- step_opt.push_back(step_opt.back());
// double minError=1000;
// double minCost=0;
// double minGamma=0;
@@ -639,10 +670,10 @@ int main(int argc, char *argv[])
double progress=0;
if(!verbose_opt[0])
pfnProgress(progress,pszMessage,pProgressArg);
- double ncost=log(ccost_opt[1])/log(step_opt[0])-log(ccost_opt[0])/log(step_opt[0]);
- double ngamma=log(gamma_opt[1])/log(step_opt[1])-log(gamma_opt[0])/log(step_opt[1]);
- for(double ccost=ccost_opt[0];ccost<=ccost_opt[1];ccost*=step_opt[0]){
- for(double gamma=gamma_opt[0];gamma<=gamma_opt[1];gamma*=step_opt[1]){
+ double ncost=log(ccost_opt[1])/log(stepcc_opt[0])-log(ccost_opt[0])/log(stepcc_opt[0]);
+ double ngamma=log(gamma_opt[1])/log(stepg_opt[0])-log(gamma_opt[0])/log(stepg_opt[0]);
+ for(double ccost=ccost_opt[0];ccost<=ccost_opt[1];ccost*=stepcc_opt[0]){
+ for(double gamma=gamma_opt[0];gamma<=gamma_opt[1];gamma*=stepg_opt[0]){
x[0]=ccost;
x[1]=gamma;
std::vector<double> theGrad;
diff --git a/src/apps/pkpolygonize.cc b/src/apps/pkpolygonize.cc
index 99fb6aa..a60d96c 100644
--- a/src/apps/pkpolygonize.cc
+++ b/src/apps/pkpolygonize.cc
@@ -39,21 +39,21 @@ using namespace std;
int main(int argc,char **argv) {
Optionpk<string> input_opt("i", "input", "Input image file");
Optionpk<string> mask_opt("m", "mask", "All pixels in the mask band with a value other than zero will be considered suitable for collection as polygons. Use input file as mask to remove background polygon! ");
- Optionpk<double> nodata_opt("nodata", "nodata", "Disgard this nodata value when creating polygons.");
Optionpk<string> output_opt("o", "output", "Output vector file");
Optionpk<string> ogrformat_opt("f", "f", "Output OGR file format","SQLite");
Optionpk<int> band_opt("b", "band", "the band to be used from input file", 0);
Optionpk<string> fname_opt("n", "name", "the field name of the output layer", "DN");
- Optionpk<short> verbose_opt("v", "verbose", "verbose mode if > 0", 0);
+ Optionpk<double> nodata_opt("nodata", "nodata", "Disgard this nodata value when creating polygons.");
+ Optionpk<short> verbose_opt("v", "verbose", "verbose mode if > 0", 0,2);
bool doProcess;//stop process when program was invoked with help option (-h --help)
try{
doProcess=input_opt.retrieveOption(argc,argv);
mask_opt.retrieveOption(argc,argv);
- nodata_opt.retrieveOption(argc,argv);
output_opt.retrieveOption(argc,argv);
ogrformat_opt.retrieveOption(argc,argv);
band_opt.retrieveOption(argc,argv);
+ nodata_opt.retrieveOption(argc,argv);
fname_opt.retrieveOption(argc,argv);
verbose_opt.retrieveOption(argc,argv);
}
@@ -62,6 +62,9 @@ int main(int argc,char **argv) {
exit(0);
}
if(!doProcess){
+ cout << endl;
+ cout << "Usage: pkpolygonize -i input [-m mask] -o output" << 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
}
diff --git a/src/apps/pkregann.cc b/src/apps/pkregann.cc
index 96a6d40..b350077 100644
--- a/src/apps/pkregann.cc
+++ b/src/apps/pkregann.cc
@@ -43,21 +43,27 @@ int main(int argc, char *argv[])
// Optionpk<float> weights_opt("w", "weights", "weights for neural network. Apply to fully connected network only, starting from first input neuron to last output neuron, including the bias neurons (last neuron in each but last layer)", 0.0);
Optionpk<float> learning_opt("l", "learning", "learning rate (default: 0.7)", 0.7);
Optionpk<unsigned int> maxit_opt("\0", "maxit", "number of maximum iterations (epoch) (default: 500)", 500);
- Optionpk<short> verbose_opt("v", "verbose", "set to: 0 (results only), 1 (confusion matrix), 2 (debug)",0);
+ Optionpk<short> verbose_opt("v", "verbose", "set to: 0 (results only), 1 (confusion matrix), 2 (debug)",0,2);
+
+ offset_opt.setHide(1);
+ scale_opt.setHide(1);
+ connection_opt.setHide(1);
+ learning_opt.setHide(1);
+ maxit_opt.setHide(1);
bool doProcess;//stop process when program was invoked with help option (-h --help)
try{
doProcess=input_opt.retrieveOption(argc,argv);
- output_opt.retrieveOption(argc,argv);
+ training_opt.retrieveOption(argc,argv);
inputCols_opt.retrieveOption(argc,argv);
outputCols_opt.retrieveOption(argc,argv);
- training_opt.retrieveOption(argc,argv);
+ output_opt.retrieveOption(argc,argv);
from_opt.retrieveOption(argc,argv);
to_opt.retrieveOption(argc,argv);
- offset_opt.retrieveOption(argc,argv);
- scale_opt.retrieveOption(argc,argv);
cv_opt.retrieveOption(argc,argv);
nneuron_opt.retrieveOption(argc,argv);
+ offset_opt.retrieveOption(argc,argv);
+ scale_opt.retrieveOption(argc,argv);
connection_opt.retrieveOption(argc,argv);
// weights_opt.retrieveOption(argc,argv);
learning_opt.retrieveOption(argc,argv);
@@ -69,6 +75,9 @@ int main(int argc, char *argv[])
exit(0);
}
if(!doProcess){
+ cout << endl;
+ cout << "Usage: pkregann -i input -t training [-ic col]* [-oc col]* -o output" << 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
}
diff --git a/src/apps/pksetmask.cc b/src/apps/pksetmask.cc
index cadf101..cc6f5a2 100644
--- a/src/apps/pksetmask.cc
+++ b/src/apps/pksetmask.cc
@@ -38,19 +38,24 @@ int main(int argc, char *argv[])
Optionpk<char> operator_opt("p", "operator", "Operator: < = > !. Use operator for each msknodata option", '=');
Optionpk<int> nodata_opt("nodata", "nodata", "nodata value to put in image if not valid", 0);
Optionpk<string> colorTable_opt("ct", "ct", "color table (file with 5 columns: id R G B ALFA (0: transparent, 255: solid)");
- Optionpk<short> verbose_opt("v", "verbose", "verbose", 0);
+ Optionpk<short> verbose_opt("v", "verbose", "verbose", 0,2);
+ otype_opt.setHide(1);
+ oformat_opt.setHide(1);
+ option_opt.setHide(1);
+ colorTable_opt.setHide(1);
+
bool doProcess;//stop process when program was invoked with help option (-h --help)
try{
doProcess=input_opt.retrieveOption(argc,argv);
mask_opt.retrieveOption(argc,argv);
+ msknodata_opt.retrieveOption(argc,argv);
output_opt.retrieveOption(argc,argv);
+ nodata_opt.retrieveOption(argc,argv);
+ operator_opt.retrieveOption(argc,argv);
otype_opt.retrieveOption(argc,argv);
oformat_opt.retrieveOption(argc,argv);
option_opt.retrieveOption(argc,argv);
- msknodata_opt.retrieveOption(argc,argv);
- operator_opt.retrieveOption(argc,argv);
- nodata_opt.retrieveOption(argc,argv);
colorTable_opt.retrieveOption(argc,argv);
verbose_opt.retrieveOption(argc,argv);
}
@@ -59,6 +64,9 @@ int main(int argc, char *argv[])
exit(0);
}
if(!doProcess){
+ cout << endl;
+ cout << "Usage: pksetmask -i input -m mask [-msknodata value] -o output" << 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
}
diff --git a/src/apps/pksieve.cc b/src/apps/pksieve.cc
index 60e3588..6dbf0de 100644
--- a/src/apps/pksieve.cc
+++ b/src/apps/pksieve.cc
@@ -42,16 +42,16 @@ int main(int argc,char **argv) {
Optionpk<string> otype_opt("ot", "otype", "Data type for output image ({Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/CInt16/CInt32/CFloat32/CFloat64}). Empty string: inherit type from input image", "");
Optionpk<string> option_opt("co", "co", "Creation option for output file. Multiple options can be specified.");
Optionpk<string> colorTable_opt("ct", "ct", "color table (file with 5 columns: id R G B ALFA (0: transparent, 255: solid)");
- Optionpk<short> verbose_opt("v", "verbose", "verbose mode if > 0", 0);
+ Optionpk<short> verbose_opt("v", "verbose", "verbose mode if > 0", 0,2);
bool doProcess;//stop process when program was invoked with help option (-h --help)
try{
doProcess=input_opt.retrieveOption(argc,argv);
- mask_opt.retrieveOption(argc,argv);
+ size_opt.retrieveOption(argc,argv);
output_opt.retrieveOption(argc,argv);
- band_opt.retrieveOption(argc,argv);
connect_opt.retrieveOption(argc,argv);
- size_opt.retrieveOption(argc,argv);
+ band_opt.retrieveOption(argc,argv);
+ mask_opt.retrieveOption(argc,argv);
otype_opt.retrieveOption(argc,argv);
option_opt.retrieveOption(argc,argv);
colorTable_opt.retrieveOption(argc,argv);
@@ -62,6 +62,9 @@ int main(int argc,char **argv) {
exit(0);
}
if(!doProcess){
+ cout << endl;
+ cout << "Usage: pksieve -i input [-s size] -o output" << 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
}
diff --git a/src/apps/pkstatascii.cc b/src/apps/pkstatascii.cc
index 1c38460..03c8fb9 100644
--- a/src/apps/pkstatascii.cc
+++ b/src/apps/pkstatascii.cc
@@ -61,17 +61,22 @@ int main(int argc, char *argv[])
Optionpk<bool> rmse_opt("rmse","rmse","calculate root mean square error between two columns (defined by -c <col1> -c <col2>",false);
Optionpk<bool> reg_opt("reg","regression","calculate linear regression between two columns and get correlation coefficient (defined by -c <col1> -c <col2>",false);
Optionpk<bool> regerr_opt("regerr","regerr","calculate linear regression between two columns and get root mean square error (defined by -c <col1> -c <col2>",false);
- Optionpk<short> verbose_opt("v", "verbose", "verbose mode when positive", 0);
+ Optionpk<short> verbose_opt("v", "verbose", "verbose mode when positive", 0,2);
+
+ src_min_opt.setHide(1);
+ src_max_opt.setHide(1);
+ fs_opt.setHide(1);
+ range_opt.setHide(1);
+ output_opt.setHide(1);
+ transpose_opt.setHide(1);
+ comment_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);
- fs_opt.retrieveOption(argc,argv);
- comment_opt.retrieveOption(argc,argv);
- output_opt.retrieveOption(argc,argv);
- transpose_opt.retrieveOption(argc,argv);
col_opt.retrieveOption(argc,argv);
- range_opt.retrieveOption(argc,argv);
+ //optional options
size_opt.retrieveOption(argc,argv);
rand_opt.retrieveOption(argc,argv);
randdist_opt.retrieveOption(argc,argv);
@@ -87,17 +92,23 @@ int main(int argc, char *argv[])
minmax_opt.retrieveOption(argc,argv);
min_opt.retrieveOption(argc,argv);
max_opt.retrieveOption(argc,argv);
- src_min_opt.retrieveOption(argc,argv);
- src_max_opt.retrieveOption(argc,argv);
histogram_opt.retrieveOption(argc,argv);
- histogram2d_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);
+ //advanced options
+ src_min_opt.retrieveOption(argc,argv);
+ src_max_opt.retrieveOption(argc,argv);
+ fs_opt.retrieveOption(argc,argv);
+ range_opt.retrieveOption(argc,argv);
+ output_opt.retrieveOption(argc,argv);
+ transpose_opt.retrieveOption(argc,argv);
+ comment_opt.retrieveOption(argc,argv);
verbose_opt.retrieveOption(argc,argv);
}
catch(string predefinedString){
@@ -105,6 +116,9 @@ int main(int argc, char *argv[])
exit(0);
}
if(!doProcess){
+ cout << endl;
+ cout << "Usage: pkstatascii -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
}
diff --git a/src/apps/pkstatogr.cc b/src/apps/pkstatogr.cc
index 88b2a1f..a31925e 100644
--- a/src/apps/pkstatogr.cc
+++ b/src/apps/pkstatogr.cc
@@ -46,13 +46,13 @@ int main(int argc, char *argv[])
Optionpk<unsigned int> nbin_opt("nbin", "nbin", "Number of bins");
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<short> verbose_opt("v", "verbose", "Verbose level", 0);
+ Optionpk<short> verbose_opt("v", "verbose", "Verbose level", 0,2);
bool doProcess;//stop process when program was invoked with help option (-h --help)
try{
doProcess=input_opt.retrieveOption(argc,argv);
- layer_opt.retrieveOption(argc,argv);
fieldname_opt.retrieveOption(argc,argv);
+ layer_opt.retrieveOption(argc,argv);
nodata_opt.retrieveOption(argc,argv);
src_min_opt.retrieveOption(argc,argv);
src_max_opt.retrieveOption(argc,argv);
@@ -74,6 +74,9 @@ int main(int argc, char *argv[])
exit(0);
}
if(!doProcess){
+ cout << endl;
+ cout << "Usage: pkstatogr -i input [-n attribute]*" << endl;
+ cout << endl;
cout << "short option -h shows basic options only, use long option --help to show all options" << endl;
exit(0);//help was invoked, stop processing
}
diff --git a/src/apps/pksvm.cc b/src/apps/pksvm.cc
index 60f9ddc..edac81a 100644
--- a/src/apps/pksvm.cc
+++ b/src/apps/pksvm.cc
@@ -94,15 +94,44 @@ int main(int argc, char *argv[])
Optionpk<unsigned int> nactive_opt("na", "nactive", "Number of active training points",1);
Optionpk<string> classname_opt("c", "class", "List of class names.");
Optionpk<short> classvalue_opt("r", "reclass", "List of class values (use same order as in class opt).");
- Optionpk<short> verbose_opt("v", "verbose", "Verbose level",0);
+ Optionpk<short> verbose_opt("v", "verbose", "Verbose level",0,2);
+
+ band_opt.setHide(1);
+ bstart_opt.setHide(1);
+ bend_opt.setHide(1);
+ balance_opt.setHide(1);
+ minSize_opt.setHide(1);
+ bag_opt.setHide(1);
+ bagSize_opt.setHide(1);
+ comb_opt.setHide(1);
+ classBag_opt.setHide(1);
+ prob_opt.setHide(1);
+ priorimg_opt.setHide(1);
+ offset_opt.setHide(1);
+ scale_opt.setHide(1);
+ svm_type_opt.setHide(1);
+ kernel_type_opt.setHide(1);
+ kernel_degree_opt.setHide(1);
+ coef0_opt.setHide(1);
+ nu_opt.setHide(1);
+ epsilon_loss_opt.setHide(1);
+ cache_opt.setHide(1);
+ epsilon_tol_opt.setHide(1);
+ shrinking_opt.setHide(1);
+ prob_est_opt.setHide(1);
+ entropy_opt.setHide(1);
+ active_opt.setHide(1);
+ nactive_opt.setHide(1);
+ verbose_opt.setHide(1);
+ random_opt.setHide(1);
bool doProcess;//stop process when program was invoked with help option (-h --help)
try{
doProcess=training_opt.retrieveOption(argc,argv);
- tlayer_opt.retrieveOption(argc,argv);
input_opt.retrieveOption(argc,argv);
output_opt.retrieveOption(argc,argv);
cv_opt.retrieveOption(argc,argv);
+ tlayer_opt.retrieveOption(argc,argv);
classname_opt.retrieveOption(argc,argv);
classvalue_opt.retrieveOption(argc,argv);
oformat_opt.retrieveOption(argc,argv);
@@ -116,6 +145,7 @@ int main(int argc, char *argv[])
mask_opt.retrieveOption(argc,argv);
msknodata_opt.retrieveOption(argc,argv);
nodata_opt.retrieveOption(argc,argv);
+ // Advanced options
band_opt.retrieveOption(argc,argv);
bstart_opt.retrieveOption(argc,argv);
bend_opt.retrieveOption(argc,argv);
@@ -150,6 +180,9 @@ int main(int argc, char *argv[])
exit(0);
}
if(!doProcess){
+ cout << endl;
+ cout << "Usage: pksvm -t training [-i input -o output] [-cv value]" << 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
}
@@ -286,11 +319,6 @@ int main(int argc, char *argv[])
cerr << "error catched" << std::endl;
exit(1);
}
- //todo: delete class 0 ?
- // if(verbose_opt[0]>=1)
- // std::cout << "erasing class 0 from training set (" << trainingMap[0].size() << " from " << totalSamples << ") samples" << std::endl;
- // totalSamples-=trainingMap[0].size();
- // trainingMap.erase(0);
//convert map to vector
// short iclass=0;
diff --git a/src/base/Optionpk.h b/src/base/Optionpk.h
index 1905978..23fdcbe 100644
--- a/src/base/Optionpk.h
+++ b/src/base/Optionpk.h
@@ -87,7 +87,8 @@ Several command line option formats are supported:
- `--longOption` (no value for boolean options, which are automatically set by invoking the option)
Option names should have regular characters and no white space in them. Some names are reserved and can not be used either:
-- short option `h` or long option `help`: shows help info
+- short option `h` or long option `help`: shows usage
+- long option `help` shows long help info
- long option `license`: shows license info
- long option `version`: shows current version of pktools
- long option `doxygen`: shows help info in table format, ready to be included in doxygen
@@ -115,6 +116,7 @@ public:
~Optionpk();
///set help information
void setHelp(const std::string& helpInfo){m_help=helpInfo;};
+ void setHide(short hide){m_hide=hide;};
bool retrieveOption(int argc, char ** argv);
template<class T1> friend std::ostream& operator<<(std::ostream & os, const Optionpk<T1>& theOption);
@@ -194,7 +196,7 @@ shortName is option invoked with `-`\n
longName is option invoked with `--`\n
helpInfo is the help message that is shown when option -h or --help is invoked\n
defaultValue is default value of the option (first value of vector: option[0])\n
-hide=0 : option is visible for in both short (`-h`) and long (`--help`) help. Typical use: mandatory options\n
+hide=0 : option is visible for in both short (`-h`). Typical use: mandatory options\n
hide=1 : option is only visible in long help (`--help`). Typical use: expert options\n
hide=2 : option is hidden for user. Typical use: Easter eggs or options only known to author
**/
@@ -356,6 +358,26 @@ template<> inline std::string string2type(std::string const& s){
return s;
}
+template<> inline double string2type(std::string const& s){
+ std::istringstream i;
+ i.precision(12);
+ i.str(s);
+ double x;
+ if (!(i >> std::setprecision(12) >> x) )
+ throw BadConversion(s);
+ return x;
+}
+
+template<> inline float string2type(std::string const& s){
+ std::istringstream i;
+ i.precision(12);
+ i.str(s);
+ float x;
+ if (!(i >> std::setprecision(12) >> x) )
+ throw BadConversion(s);
+ return x;
+}
+
///specialization for OGRFieldType
template<> inline OGRFieldType string2type(std::string const& s){
OGRFieldType ftype;
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/pkg-grass/pktools.git
More information about the Pkg-grass-devel
mailing list