[pktools] 280/375: added pkkalman.cc
Bas Couwenberg
sebastic at xs4all.nl
Wed Dec 3 21:54:22 UTC 2014
This is an automated email from the git hooks/post-receive script.
sebastic-guest pushed a commit to branch upstream-master
in repository pktools.
commit 09265f1663f44de48da33939e1d84b91ac75d420
Author: user <user at osgeolive.(none)>
Date: Tue Jun 10 13:03:32 2014 +0200
added pkkalman.cc
---
src/algorithms/StatFactory.h | 14 +++-
src/apps/Makefile.am | 4 +-
src/apps/pkkalman.cc | 194 +++++++++++++++++++++++++++++++++++++++++++
src/apps/pkstatascii.cc | 16 +++-
4 files changed, 224 insertions(+), 4 deletions(-)
diff --git a/src/algorithms/StatFactory.h b/src/algorithms/StatFactory.h
index e6c410b..d8f76a7 100644
--- a/src/algorithms/StatFactory.h
+++ b/src/algorithms/StatFactory.h
@@ -178,6 +178,7 @@ public:
template<class T> double correlation(const std::vector<T>& x, const std::vector<T>& y, int delay=0) const;
template<class T> double cross_correlation(const std::vector<T>& x, const std::vector<T>& y, int maxdelay, std::vector<T>& z) const;
template<class T> double linear_regression(const std::vector<T>& x, const std::vector<T>& y, double &c0, double &c1) const;
+ template<class T> double linear_regression_err(const std::vector<T>& x, const std::vector<T>& y, double &c0, double &c1) const;
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);
@@ -871,7 +872,7 @@ template<class T> double StatFactory::cross_correlation(const std::vector<T>& x,
return sumCorrelation;
}
- template<class T> double StatFactory::linear_regression(const std::vector<T>& x, const std::vector<T>& y, double &c0, double &c1) const{
+template<class T> double StatFactory::linear_regression(const std::vector<T>& x, const std::vector<T>& y, double &c0, double &c1) const{
assert(x.size()==y.size());
assert(x.size());
double cov00;
@@ -882,6 +883,17 @@ template<class T> double StatFactory::cross_correlation(const std::vector<T>& x,
return (1-sumsq/var(y)/(y.size()-1));
}
+template<class T> double StatFactory::linear_regression_err(const std::vector<T>& x, const std::vector<T>& y, double &c0, double &c1) const{
+ assert(x.size()==y.size());
+ assert(x.size());
+ double cov00;
+ double cov01;
+ double cov11;
+ double sumsq;
+ gsl_fit_linear(&(x[0]),1,&(y[0]),1,x.size(),&c0,&c1,&cov00,&cov01,&cov11,&sumsq);
+ return sqrt((sumsq)/(y.size()));
+}
+
//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)
diff --git a/src/apps/Makefile.am b/src/apps/Makefile.am
index b2d67ac..7fd00be 100644
--- a/src/apps/Makefile.am
+++ b/src/apps/Makefile.am
@@ -6,7 +6,7 @@ LDADD = $(GSL_LIBS) $(GDAL_LDFLAGS) $(top_builddir)/src/algorithms/libalgorithms
###############################################################################
# the program to build and install (the names of the final binaries)
-bin_PROGRAMS = pkinfo pkcrop pkreclass pkgetmask pksetmask pkcreatect pkdumpimg pkdumpogr pksieve pkstatascii pkstatogr pkegcs pkextract pkfillnodata pkfilter pkfilterdem pkenhance pkfilterascii pkdsm2shadow pkcomposite pkndvi pkpolygonize pkascii2img pksvm pkfssvm pkascii2ogr pkeditogr
+bin_PROGRAMS = pkinfo pkcrop pkreclass pkgetmask pksetmask pkcreatect pkdumpimg pkdumpogr pksieve pkstatascii pkstatogr pkegcs pkextract pkfillnodata pkfilter pkkalman pkfilterdem pkenhance pkfilterascii pkdsm2shadow pkcomposite pkndvi pkpolygonize pkascii2img pksvm pkfssvm pkascii2ogr pkeditogr
# the program to build but not install (the names of the final binaries)
#noinst_PROGRAMS = pkxcorimg pkgeom
@@ -53,6 +53,8 @@ pkextract_SOURCES = pkextract.cc
pkfillnodata_SOURCES = pkfillnodata.cc
pkfilter_SOURCES = pkfilter.cc
pkfilter_LDADD = $(GSL_LIBS) $(AM_LDFLAGS) -lgsl -lgdal
+pkkalman_SOURCES = pkkalman.cc
+pkkalman_LDADD = $(GSL_LIBS) $(AM_LDFLAGS) -lgsl -lgdal
pkfilterdem_SOURCES = pkfilterdem.cc
pkenhance_SOURCES = pkenhance.cc
pkenhance_LDADD = $(AM_LDFLAGS) -lgdal
diff --git a/src/apps/pkkalman.cc b/src/apps/pkkalman.cc
new file mode 100644
index 0000000..df9abf1
--- /dev/null
+++ b/src/apps/pkkalman.cc
@@ -0,0 +1,194 @@
+/**********************************************************************
+pkkalman.cc: program to kalman raster images: median, min/max, morphological, kalmaning
+Copyright (C) 2008-2014 Pieter Kempeneers
+
+This file is part of pktools
+
+pktools is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+pktools is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with pktools. If not, see <http://www.gnu.org/licenses/>.
+***********************************************************************/
+#include <sstream>
+#include <vector>
+#include "base/Optionpk.h"
+#include "base/Vector2d.h"
+#include "algorithms/StatFactory.h"
+#include "imageclasses/ImgReaderGdal.h"
+#include "imageclasses/ImgWriterGdal.h"
+
+using namespace std;
+/*------------------
+ Main procedure
+ ----------------*/
+int main(int argc,char **argv) {
+ Optionpk<string> model_opt("mod","model","model input datasets, e.g., MODIS (use: -mod model1 -mod model2 etc.");
+ Optionpk<string> observation_opt("obs","observation","observation input datasets, e.g., landsat (use: -obs obs1 -obs obs2 etc.");
+ Optionpk<int> tobservation_opt("tobs","tobservation","time sequence of observation input (sequence must have exact same length as observation input)");
+ Optionpk<string> output_opt("o", "output", "Suffix for output image datasets");
+ Optionpk<float> threshold_opt("t", "threshold", "threshold for selecting samples (randomly). Provide probability in percentage (>0) or absolute (<0).", 100);
+ Optionpk<double> modnodata_opt("modnodata", "modnodata", "invalid value for model input", 0);
+ Optionpk<int> bndmodnodata_opt("bmnodata", "bndmodnodata", "Bands in model input to check if pixel is valid (used for srcnodata, min and max options)", 0);
+ Optionpk<double> obsnodata_opt("obsnodata", "obsnodata", "invalid value for observation input", 0);
+ Optionpk<int> bndobsnodata_opt("bonodata", "bndobsnodata", "Bands in observation input to check if pixel is valid (used for srcnodata, min and max options)", 0);
+ Optionpk<int> down_opt("down", "down", "Downsampling factor for reading model data to calculate regression", 90);
+ Optionpk<short> verbose_opt("v", "verbose", "verbose mode when positive", 0);
+
+ bool doProcess;//stop process when program was invoked with help option (-h --help)
+ try{
+ doProcess=model_opt.retrieveOption(argc,argv);
+ observation_opt.retrieveOption(argc,argv);
+ tobservation_opt.retrieveOption(argc,argv);
+ output_opt.retrieveOption(argc,argv);
+ threshold_opt.retrieveOption(argc,argv);
+ modnodata_opt.retrieveOption(argc,argv);
+ bndmodnodata_opt.retrieveOption(argc,argv);
+ obsnodata_opt.retrieveOption(argc,argv);
+ bndobsnodata_opt.retrieveOption(argc,argv);
+ down_opt.retrieveOption(argc,argv);
+ verbose_opt.retrieveOption(argc,argv);
+ }
+ catch(string predefinedString){
+ std::cout << predefinedString << std::endl;
+ exit(0);
+ }
+ if(!doProcess){
+ 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
+ }
+
+ try{
+ ostringstream errorStream;
+ if(model_opt.size()<2){
+ errorStream << "Error: no model dataset selected, use option -mod" << endl;
+ throw(errorStream.str());
+ }
+ if(observation_opt.size()<1){
+ errorStream << "Error: no observation dataset selected, use option -obs" << endl;
+ throw(errorStream.str());
+ }
+ if(model_opt.size()<observation_opt.size()){
+ errorStream << "Error: sequence of models should be larger than observations" << endl;
+ throw(errorStream.str());
+ }
+ if(tobservation_opt.size()!=observation_opt.size()){
+ errorStream << "Error: time sequence for observation must match size of observation dataset" << endl;
+ throw(errorStream.str());
+ }
+ if(output_opt.empty()){
+ errorStream << "Error: suffix for output datasets is not provided" << endl;
+ throw(errorStream.str());
+ }
+ }
+ catch(string errorString){
+ std::cout << errorString << std::endl;
+ exit(1);
+ }
+
+ while(modnodata_opt.size()<bndmodnodata_opt.size())
+ modnodata_opt.push_back(modnodata_opt[0]);
+ while(bndmodnodata_opt.size()<modnodata_opt.size())
+ bndmodnodata_opt.push_back(bndmodnodata_opt[0]);
+ while(obsnodata_opt.size()<bndobsnodata_opt.size())
+ obsnodata_opt.push_back(obsnodata_opt[0]);
+ while(bndobsnodata_opt.size()<obsnodata_opt.size())
+ bndobsnodata_opt.push_back(bndobsnodata_opt[0]);
+
+ statfactory::StatFactory stat;
+
+ vector<ImgReaderGdal> imgReaderModel(model_opt.size());
+ vector<ImgReaderGdal> imgReaderObs(observation_opt.size());
+ vector<ImgWriterGdal> imgWriterPred(model_opt.size());
+
+ int obsindex=0;
+ //forward step
+
+ const char* pszMessage;
+ void* pProgressArg=NULL;
+ GDALProgressFunc pfnProgress=GDALTermProgress;
+ double progress=0;
+ srand(time(NULL));
+
+ for(int modindex=0;modindex<model_opt.size()-1;++modindex){
+ //calculate regression between two subsequence model inputs
+ imgReaderModel[modindex].open(model_opt[modindex]);
+ imgReaderModel[modindex+1].open(model_opt[modindex]);
+ //calculate regression
+ //we could re-use the points from second image from last run, but
+ //to keep it general, we must redo it (overlap might have changed)
+ pfnProgress(progress,pszMessage,pProgressArg);
+ int icol1=0,irow1=0;
+ vector<double> rowBuffer1(imgReaderModel[modindex].nrOfCol());
+ vector<double> rowBuffer2(imgReaderModel[modindex+1].nrOfCol());
+ vector<double> buffer1;
+ vector<double> buffer2;
+
+ for(irow1=0;irow1<imgReaderModel[modindex].nrOfRow();++irow1){
+ if(irow1%down_opt[0])
+ continue;
+ icol1=0;
+ double icol2=0,irow2=0;
+ double geox=0,geoy=0;
+ imgReaderModel[modindex].readData(rowBuffer1,GDT_Float64,irow1);
+ imgReaderModel[modindex].image2geo(icol1,irow1,geox,geoy);
+ imgReaderModel[modindex+1].geo2image(geox,geoy,icol2,irow2);
+ icol2=static_cast<int>(icol2);
+ irow2=static_cast<int>(irow2);
+ imgReaderModel[modindex+1].readData(rowBuffer2,GDT_Float64,irow2);
+ for(icol1=0;icol1<imgReaderModel[modindex].nrOfCol();++icol1){
+ if(icol1%down_opt[0])
+ continue;
+ if(threshold_opt[0]>0){//percentual value
+ double p=static_cast<double>(rand())/(RAND_MAX);
+ p*=100.0;
+ if(p>threshold_opt[0])
+ continue;//do not select for now, go to next column
+ }
+ else if(buffer1.size()>-threshold_opt[0])//absolute value
+ continue;//do not select any more pixels
+ imgReaderModel[modindex].image2geo(icol1,irow1,geox,geoy);
+ imgReaderModel[modindex+1].geo2image(geox,geoy,icol2,irow2);
+ icol2=static_cast<int>(icol2);
+ irow2=static_cast<int>(irow2);
+ //check for nodata
+ double valmod1=rowBuffer1[icol1];
+ double valmod2=rowBuffer2[icol2];
+ bool readValid=true;
+ for(int vband=0;vband<bndmodnodata_opt.size();++vband){
+ if(modnodata_opt.size()>vband){
+ if(valmod1==modnodata_opt[vband] || valmod2==modnodata_opt[vband]){
+ readValid=false;
+ break;
+ }
+ }
+ }
+ buffer1.push_back(valmod1);
+ buffer2.push_back(valmod2);
+ }
+ }
+ double c0=0;
+ double c1=0;
+ double err=stat.linear_regression_err(buffer1,buffer2,c0,c1);
+ if(verbose_opt[0])
+ cout << "linear regression model-model based on " << buffer1.size() << " points: " << c0 << "+" << c1 << " * x " << " with rmse: " << err << endl;
+
+ if(tobservation_opt[obsindex]==modindex){//update
+ imgReaderObs[obsindex].open(observation_opt[obsindex]);
+ //calculate regression
+ //set ImgWriterPred[modindex]=imgReaderObs[obindex]
+ //calculate regression between model and observation
+
+ }
+ else{//prediction
+ }
+ ++obsindex;
+ }
+}
diff --git a/src/apps/pkstatascii.cc b/src/apps/pkstatascii.cc
index 69a1137..fe3b842 100644
--- a/src/apps/pkstatascii.cc
+++ b/src/apps/pkstatascii.cc
@@ -59,8 +59,9 @@ int main(int argc, char *argv[])
Optionpk<bool> kde_opt("kde","kde","Use Kernel density estimation when producing histogram. The standard deviation is estimated based on Silverman's rule of thumb",false);
Optionpk<bool> correlation_opt("cor","correlation","calculate Pearson produc-moment correlation coefficient between two columns (defined by -c <col1> -c <col2>",false);
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 error between two columns (defined by -c <col1> -c <col2>",false);
- Optionpk<short> verbose_opt("v", "verbose", "verbose mode when > 0", 0);
+ 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);
bool doProcess;//stop process when program was invoked with help option (-h --help)
try{
@@ -96,6 +97,7 @@ int main(int argc, char *argv[])
correlation_opt.retrieveOption(argc,argv);
rmse_opt.retrieveOption(argc,argv);
reg_opt.retrieveOption(argc,argv);
+ regerr_opt.retrieveOption(argc,argv);
verbose_opt.retrieveOption(argc,argv);
}
catch(string predefinedString){
@@ -292,6 +294,16 @@ int main(int argc, char *argv[])
double r2=stat.linear_regression(dataVector[0],dataVector[1],c0,c1);
cout << "linear regression between columns: " << col_opt[0] << " and " << col_opt[1] << ": " << c0 << "+" << c1 << " * x " << " with R^2 (square correlation coefficient): " << r2 << endl;
}
+ if(regerr_opt[0]){
+ assert(dataVector.size()==2);
+ double c0=0;
+ double c1=0;
+ double err=stat.linear_regression_err(dataVector[0],dataVector[1],c0,c1);
+ if(verbose_opt[0])
+ cout << "linear regression between columns: " << col_opt[0] << " and " << col_opt[1] << ": " << c0 << "+" << c1 << " * x " << " with rmse: " << err << endl;
+ else
+ cout << c0 << " " << c1 << " " << err << endl;
+ }
if(histogram_opt[0]){
for(int irow=0;irow<statVector.begin()->size();++irow){
double binValue=0;
--
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