[pktools] 281/375: added class ImgRegression
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 ea06516d71792a8228924848662be1e41ebbc45b
Author: user <user at osgeolive.(none)>
Date: Tue Jun 10 15:19:34 2014 +0200
added class ImgRegression
---
src/algorithms/ImgRegression.cc | 82 ++++++++++++++++++++++++++
src/algorithms/ImgRegression.h | 42 ++++++++++++++
src/apps/Makefile.am | 4 +-
src/apps/pkkalman.cc | 121 ++++++++++++++++++---------------------
src/imageclasses/ImgReaderGdal.h | 1 +
5 files changed, 184 insertions(+), 66 deletions(-)
diff --git a/src/algorithms/ImgRegression.cc b/src/algorithms/ImgRegression.cc
new file mode 100644
index 0000000..3ebf557
--- /dev/null
+++ b/src/algorithms/ImgRegression.cc
@@ -0,0 +1,82 @@
+/**********************************************************************
+ImgRegression.cc: class to calculate regression between two raster datasets
+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 "ImgRegression.h"
+#include <iostream>
+
+using namespace imgregression;
+
+ImgRegression::ImgRegression(void)
+: m_threshold(0), m_down(1)
+{}
+
+ImgRegression::~ImgRegression(void)
+{}
+
+double ImgRegression::getRMSE(const ImgReaderGdal& imgReader1, const ImgReaderGdal& imgReader2, double& c0, double& c1, short verbose) const{
+ int icol1=0,irow1=0;
+ std::vector<double> rowBuffer1(imgReader1.nrOfCol());
+ std::vector<double> rowBuffer2(imgReader2.nrOfCol());
+ std::vector<double> buffer1;
+ std::vector<double> buffer2;
+
+ for(irow1=0;irow1<imgReader1.nrOfRow();++irow1){
+ if(irow1%m_down)
+ continue;
+ icol1=0;
+ double icol2=0,irow2=0;
+ double geox=0,geoy=0;
+ imgReader1.readData(rowBuffer1,GDT_Float64,irow1);
+ imgReader1.image2geo(icol1,irow1,geox,geoy);
+ imgReader2.geo2image(geox,geoy,icol2,irow2);
+ icol2=static_cast<int>(icol2);
+ irow2=static_cast<int>(irow2);
+ imgReader2.readData(rowBuffer2,GDT_Float64,irow2);
+ for(icol1=0;icol1<imgReader1.nrOfCol();++icol1){
+ if(icol1%m_down)
+ continue;
+ if(m_threshold>0){//percentual value
+ double p=static_cast<double>(rand())/(RAND_MAX);
+ p*=100.0;
+ if(p>m_threshold)
+ continue;//do not select for now, go to next column
+ }
+ imgReader1.image2geo(icol1,irow1,geox,geoy);
+ imgReader2.geo2image(geox,geoy,icol2,irow2);
+ icol2=static_cast<int>(icol2);
+ irow2=static_cast<int>(irow2);
+ //check for nodata
+ double value1=rowBuffer1[icol1];
+ double value2=rowBuffer2[icol2];
+ bool readValid=true;
+ if(imgReader1.isNoData(value1)||imgReader2.isNoData(value2))
+ readValid=false;
+
+ buffer1.push_back(value1);
+ buffer2.push_back(value2);
+ if(verbose>1)
+ std::cout << geox << " " << geoy << " " << icol1 << " " << irow1 << " " << icol2 << " " << irow2 << " " << buffer1.back() << " " << buffer2.back() << std::endl;
+ }
+ }
+ statfactory::StatFactory stat;
+ double err=stat.linear_regression_err(buffer1,buffer2,c0,c1);
+ if(verbose)
+ std::cout << "linear regression model-model based on " << buffer1.size() << " points: " << c0 << "+" << c1 << " * x " << " with rmse: " << err << std::endl;
+ return err;
+}
diff --git a/src/algorithms/ImgRegression.h b/src/algorithms/ImgRegression.h
new file mode 100644
index 0000000..12ba62d
--- /dev/null
+++ b/src/algorithms/ImgRegression.h
@@ -0,0 +1,42 @@
+/**********************************************************************
+ImgRegression.h: class to calculate regression between two raster datasets
+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/>.
+***********************************************************************/
+#ifndef _IMGREGRESSION_H_
+#define _IMGREGRESSION_H_
+
+#include <vector>
+#include "imageclasses/ImgReaderGdal.h"
+#include "imageclasses/ImgWriterGdal.h"
+#include "StatFactory.h"
+
+namespace imgregression
+{
+ class ImgRegression{
+ public:
+ ImgRegression(void);
+ ~ImgRegression(void);
+ double getRMSE(const ImgReaderGdal& imgReader1, const ImgReaderGdal& imgReader2, double &c0, double &c1, short verbose=0) const;
+ void setThreshold(double theThreshold){m_threshold=theThreshold;};
+ void setDown(int theDown){m_down=theDown;};
+ private:
+ int m_down;
+ double m_threshold;
+ };
+}
+#endif //_IMGREGRESSION_H_
diff --git a/src/apps/Makefile.am b/src/apps/Makefile.am
index 7fd00be..1d4875e 100644
--- a/src/apps/Makefile.am
+++ b/src/apps/Makefile.am
@@ -53,8 +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
+pkkalman_SOURCES = pkkalman.cc $(top_srcdir)/src/algorithms/ImgRegression.h $(top_srcdir)/src/algorithms/ImgRegression.cc
+pkkalman_LDADD = -lalgorithms $(GSL_LIBS) $(AM_LDFLAGS) -lalgorithms
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
index df9abf1..38afbe5 100644
--- a/src/apps/pkkalman.cc
+++ b/src/apps/pkkalman.cc
@@ -21,9 +21,10 @@ along with pktools. If not, see <http://www.gnu.org/licenses/>.
#include <vector>
#include "base/Optionpk.h"
#include "base/Vector2d.h"
-#include "algorithms/StatFactory.h"
#include "imageclasses/ImgReaderGdal.h"
#include "imageclasses/ImgWriterGdal.h"
+#include "algorithms/StatFactory.h"
+#include "algorithms/ImgRegression.h"
using namespace std;
/*------------------
@@ -34,11 +35,9 @@ int main(int argc,char **argv) {
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<float> threshold_opt("t", "threshold", "threshold for selecting samples (randomly). Provide probability in percentage (>0) or absolute (<0).", 0);
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);
@@ -50,9 +49,7 @@ int main(int argc,char **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);
}
@@ -93,17 +90,7 @@ int main(int argc,char **argv) {
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;
-
+ imgregression::ImgRegression imgreg;
vector<ImgReaderGdal> imgReaderModel(model_opt.size());
vector<ImgReaderGdal> imgReaderObs(observation_opt.size());
vector<ImgWriterGdal> imgWriterPred(model_opt.size());
@@ -117,10 +104,15 @@ int main(int argc,char **argv) {
double progress=0;
srand(time(NULL));
+ imgreg.setDown(down_opt[0]);
+ imgreg.setThreshold(threshold_opt[0]);
+
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]);
+ imgReaderModel[modindex+1].open(model_opt[modindex+1]);
+ imgReaderModel[modindex].setNoData(modnodata_opt);
+ imgReaderModel[modindex+1].setNoData(modnodata_opt);
//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)
@@ -131,54 +123,54 @@ int main(int argc,char **argv) {
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;
+
+ // 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
+ // }
+ // 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);
+ // if(verbose_opt[0]>1)
+ // cout << geox << " " << geoy << " " << icol1 << " " << irow1 << " " << icol2 << " " << irow2 << " " << buffer1.back() << " " << buffer2.back() << endl;
+ // }
+ // }
+
+ double err=imgreg.getRMSE(imgReaderModel[modindex],imgReaderModel[modindex+1],c0,c1,verbose_opt[0]);
if(tobservation_opt[obsindex]==modindex){//update
imgReaderObs[obsindex].open(observation_opt[obsindex]);
@@ -192,3 +184,4 @@ int main(int argc,char **argv) {
++obsindex;
}
}
+
diff --git a/src/imageclasses/ImgReaderGdal.h b/src/imageclasses/ImgReaderGdal.h
index 554e615..5b68630 100644
--- a/src/imageclasses/ImgReaderGdal.h
+++ b/src/imageclasses/ImgReaderGdal.h
@@ -80,6 +80,7 @@ public:
int getNoDataValues(std::vector<double>& noDataValues) const;
bool isNoData(double value) const{if(m_noDataValues.empty()) return false;else return find(m_noDataValues.begin(),m_noDataValues.end(),value)!=m_noDataValues.end();};
int pushNoDataValue(double noDataValue);
+ int setNoData(const std::vector<double> nodata){m_noDataValues=nodata;};
CPLErr GDALSetNoDataValue(double noDataValue, int band=0) {return getRasterBand(band)->SetNoDataValue(noDataValue);};
bool covers(double x, double y) const;
bool covers(double ulx, double uly, double lrx, double lry) const;
--
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