[pktools] 138/375: moved pkstat.cc to pkstatascii.cc
Bas Couwenberg
sebastic at xs4all.nl
Wed Dec 3 21:54:07 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 67e2c1117f366dc078d40e7fa5f6e898c2074020
Author: Pieter Kempeneers <kempenep at gmail.com>
Date: Fri Nov 15 18:12:10 2013 +0100
moved pkstat.cc to pkstatascii.cc
---
ChangeLog | 8 +-
src/algorithms/ConfusionMatrix.cc | 8 +-
src/algorithms/ConfusionMatrix.h | 6 +-
src/algorithms/FeatureSelector.h | 40 ++--
src/algorithms/Filter.cc | 55 +++++
src/algorithms/Filter.h | 2 +-
src/algorithms/StatFactory.h | 365 +++++++++++++++++++++++----------
src/algorithms/myfann_cpp.h | 16 +-
src/apps/Makefile.am | 6 +-
src/apps/pkegcs.cc | 20 +-
src/apps/pkinfo.cc | 38 ++--
src/apps/pkregression_nn.cc | 2 +
src/apps/{pkstat.cc => pkstatascii.cc} | 100 ++++++++-
src/apps/pkstatogr.cc | 24 +--
src/base/PosValue.h | 2 -
src/base/Vector2d.h | 10 +-
src/imageclasses/ImgReaderOgr.h | 130 ++++++------
17 files changed, 557 insertions(+), 275 deletions(-)
diff --git a/ChangeLog b/ChangeLog
index baf7b33..1778328 100755
--- a/ChangeLog
+++ b/ChangeLog
@@ -117,7 +117,10 @@ version 2.4.1
todo: take priors into account for cross validation
ordering of labels before training
version 2.4.2
- - todo: change all projection options to a_srs
+ - general
+ removed using namespace std from header files
+ - PosValue.h
+ remove using namespace std;
- apps/Makefile.am
add GSL_LIBS to AM_LDFLAGS and LDADD
todo: remove redundancy in AM_LDFLAGS and LDADD
@@ -137,3 +140,6 @@ version 2.4.2
correct for wrong include path for FileReaderLas.h (lasclasses instead of fileclasses)
- pkfillnodata
default maximum distance changed from 3 to 0 (infinity)
+ - still todo
+ change all projection options to a_srs
+ rename pkstat to pkstatascii
diff --git a/src/algorithms/ConfusionMatrix.cc b/src/algorithms/ConfusionMatrix.cc
index 02a4653..64b2945 100644
--- a/src/algorithms/ConfusionMatrix.cc
+++ b/src/algorithms/ConfusionMatrix.cc
@@ -66,16 +66,16 @@ ConfusionMatrix& ConfusionMatrix::operator=(const ConfusionMatrix& cm){
ConfusionMatrix& ConfusionMatrix::operator+=(const ConfusionMatrix &cm)
{
if(cm.m_classes.size()!=this->m_classes.size()){
- cerr << "error0: "<< cm.m_classes.size() << "!=" << this->m_classes.size() << endl;
+ std::cerr << "error0: "<< cm.m_classes.size() << "!=" << this->m_classes.size() << std::endl;
exit(0);
}
if(cm.m_results.size()!=this->m_results.size()){
- cerr << "error1: "<< cm.m_results.size() << "!=" << this->m_results.size() << endl;
+ std::cerr << "error1: "<< cm.m_results.size() << "!=" << this->m_results.size() << std::endl;
exit(1);
}
for(int irow=0;irow<m_results.size();++irow){
if(cm.m_results[irow].size()!=this->m_results[irow].size()){
- cerr << "error2: " << cm.m_results[irow].size() << "!=" << this->m_results[irow].size() << endl;
+ std::cerr << "error2: " << cm.m_results[irow].size() << "!=" << this->m_results[irow].size() << std::endl;
exit(2);
}
for(int icol=0;icol<m_results[irow].size();++icol)
@@ -157,7 +157,7 @@ void ConfusionMatrix::incrementResult(const std::string& theRef, const std::stri
int ic=distance(m_classes.begin(),find(m_classes.begin(),m_classes.end(),theClass));
assert(ir>=0);
if(ir>=m_results.size())
- cerr << "Error: " << theRef << " not found in class ConfusionMatrix when incrementing for class " << theClass << endl;
+ std::cerr << "Error: " << theRef << " not found in class ConfusionMatrix when incrementing for class " << theClass << std::endl;
assert(ir<m_results.size());
assert(ic>=0);
assert(ic<m_results[ir].size());
diff --git a/src/algorithms/ConfusionMatrix.h b/src/algorithms/ConfusionMatrix.h
index fd3edc0..a6ae270 100644
--- a/src/algorithms/ConfusionMatrix.h
+++ b/src/algorithms/ConfusionMatrix.h
@@ -69,16 +69,16 @@ public:
return ConfusionMatrix(*this)+=cm;
}
void sortClassNames();
- friend ostream& operator<<(ostream& os, const ConfusionMatrix &cm){
+ friend std::ostream& operator<<(std::ostream& os, const ConfusionMatrix &cm){
for(int iclass=0;iclass<cm.nClasses();++iclass)
os << "\t" << cm.m_classes[iclass];
- os << endl;
+ os << std::endl;
assert(cm.m_classes.size()==cm.m_results.size());
for(int irow=0;irow<cm.m_results.size();++irow){
os << cm.m_classes[irow];
for(int icol=0;icol<cm.m_results[irow].size();++icol)
os << "\t" << cm.m_results[irow][icol];
- os << endl;
+ os << std::endl;
}
return os;
};
diff --git a/src/algorithms/FeatureSelector.h b/src/algorithms/FeatureSelector.h
index e40cdb2..fafd98a 100644
--- a/src/algorithms/FeatureSelector.h
+++ b/src/algorithms/FeatureSelector.h
@@ -85,7 +85,7 @@ template<class T> double FeatureSelector::forwardUnivariate(std::vector< Vector2
if(cost[ilevel].value>0)
subset.push_back(cost[ilevel].position);
if(verbose)
- cout << "feature " << subset.back() << " has cost: " << cost[ilevel].value << endl;
+ std::cout << "feature " << subset.back() << " has cost: " << cost[ilevel].value << std::endl;
++ilevel;
}
double maxCost=-1;
@@ -116,9 +116,9 @@ template<class T> double FeatureSelector::forward(std::vector< Vector2d<T> >& v,
maxCost=addFeature(v,*getCost,subset,verbose);
if(verbose){
for(std::list<int>::const_iterator lit=subset.begin();lit!=subset.end();++lit)
- cout << *lit << " ";
- cout << endl;
- // cout << "added " << subset.back() << ", " << subset.size() << "/" << maxFeatures << " features selected with cost: " << maxCost << endl;
+ std::cout << *lit << " ";
+ std::cout << std::endl;
+ // std::cout << "added " << subset.back() << ", " << subset.size() << "/" << maxFeatures << " features selected with cost: " << maxCost << std::endl;
}
}//while
return maxCost;
@@ -139,9 +139,9 @@ template<class T> double FeatureSelector::backward(std::vector< Vector2d<T> >& v
maxCost=removeFeature(v,*getCost,subset,removedFeature,verbose);
if(verbose){
for(std::list<int>::const_iterator lit=subset.begin();lit!=subset.end();++lit)
- cout << *lit << " ";
- cout << endl;
- // cout << "removed " << removedFeature << ", " << subset.size() << "/" << minFeatures << " features remain with cost: " << maxCost << endl;
+ std::cout << *lit << " ";
+ std::cout << std::endl;
+ // std::cout << "removed " << removedFeature << ", " << subset.size() << "/" << minFeatures << " features remain with cost: " << maxCost << std::endl;
}
}//while
return maxCost;
@@ -161,21 +161,21 @@ template<class T> double FeatureSelector::floating(std::vector< Vector2d<T> >& v
cost.push_back(addFeature(v,*getCost,subset,verbose));
++k;
if(verbose>1)
- cout << "added " << subset.back() << ", " << cost.size()-1 << "/" << maxFeatures << " features selected with cost: " << cost.back() << endl;
+ std::cout << "added " << subset.back() << ", " << cost.size()-1 << "/" << maxFeatures << " features selected with cost: " << cost.back() << std::endl;
else if(verbose){
for(std::list<int>::const_iterator lit=subset.begin();lit!=subset.end();++lit)
- cout << *lit << " ";
- cout << endl;
+ std::cout << *lit << " ";
+ std::cout << std::endl;
}
while(k<maxFeatures){
cost.push_back(addFeature(v,*getCost,subset,verbose));
++k;
if(verbose>1)
- cout << "added " << subset.back() << ", " << cost.size()-1 << "/" << maxFeatures << " features selected with cost: " << cost.back() << endl;
+ std::cout << "added " << subset.back() << ", " << cost.size()-1 << "/" << maxFeatures << " features selected with cost: " << cost.back() << std::endl;
else if(verbose){
for(std::list<int>::const_iterator lit=subset.begin();lit!=subset.end();++lit)
- cout << *lit << " ";
- cout << " (cost: " << cost.back() << ")" << endl;
+ std::cout << *lit << " ";
+ std::cout << " (cost: " << cost.back() << ")" << std::endl;
}
while(k>1){
@@ -186,11 +186,11 @@ template<class T> double FeatureSelector::floating(std::vector< Vector2d<T> >& v
cost[k]=cost_r;
cost.pop_back();
if(verbose>1)
- cout << "removed " << x_r << ", " << cost.size()-1 << "/" << maxFeatures << " features remain with cost: " << cost_r << endl;
+ std::cout << "removed " << x_r << ", " << cost.size()-1 << "/" << maxFeatures << " features remain with cost: " << cost_r << std::endl;
else if(verbose){
for(std::list<int>::const_iterator lit=subset.begin();lit!=subset.end();++lit)
- cout << *lit << " ";
- cout << " (cost: " << cost.back() << ")" << endl;
+ std::cout << *lit << " ";
+ std::cout << " (cost: " << cost.back() << ")" << std::endl;
}
continue;
}
@@ -199,7 +199,7 @@ template<class T> double FeatureSelector::floating(std::vector< Vector2d<T> >& v
break;
}
else if(verbose)
- cout << "could not remove any feature" << endl;
+ std::cout << "could not remove any feature" << std::endl;
cost.pop_back();
}
}
@@ -238,7 +238,7 @@ template<class T> double FeatureSelector::bruteForce(std::vector< Vector2d<T> >&
catch(...){ //singular matrix encountered
catchset=tmpset;//this tmpset resulted in failure of getCost
if(verbose){
- cout << "Could not get cost from set: " << endl;
+ std::cout << "Could not get cost from set: " << std::endl;
gsl_combination_fprintf(stdout,c," %u");
printf("\n");
}
@@ -288,7 +288,7 @@ template<class T> double FeatureSelector::addFeature(std::vector< Vector2d<T> >&
catch(...){
catchset=tmpset;//this tmpset resulted in singular matrix
if(verbose)
- cout << "Could not add feature " << tmpset.back() << endl;
+ std::cout << "Could not add feature " << tmpset.back() << std::endl;
tmpset.pop_back();
continue;
}
@@ -336,7 +336,7 @@ template<class T> double FeatureSelector::removeFeature(std::vector< Vector2d<T>
catch(...){
catchset=tmpset;//this tmpset resulted in singular matrix
if(verbose)
- cout << "Could not remove feature " << last << endl;
+ std::cout << "Could not remove feature " << last << std::endl;
tmpset.push_front(last);
continue;
}
diff --git a/src/algorithms/Filter.cc b/src/algorithms/Filter.cc
index 56d3bf0..69812c7 100644
--- a/src/algorithms/Filter.cc
+++ b/src/algorithms/Filter.cc
@@ -131,6 +131,61 @@ void filter::Filter::doit(const ImgReaderGdal& input, ImgWriterGdal& output, sho
}
}
+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
+ int nband=srf[0].size();
+ double start=floor(wavelengthIn[0]);
+ double end=ceil(wavelengthIn.back());
+ if(verbose)
+ std::cout << "wavelengths in [" << start << "," << end << "]" << std::endl << std::flush;
+
+ statfactory::StatFactory stat;
+
+ gsl_interp_accel *acc;
+ stat.allocAcc(acc);
+ gsl_spline *spline;
+ stat.getSpline(interpolationType,nband,spline);
+ stat.initSpline(spline,&(srf[0][0]),&(srf[1][0]),nband);
+ if(verbose)
+ std::cout << "calculating norm of srf" << std::endl << std::flush;
+ double norm=0;
+ norm=gsl_spline_eval_integ(spline,srf[0].front(),srf[0].back(),acc);
+ if(verbose)
+ std::cout << "norm of srf: " << norm << std::endl << std::flush;
+ gsl_spline_free(spline);
+ gsl_interp_accel_free(acc);
+ std::vector<double> wavelength_fine;
+ for(double win=floor(wavelengthIn[0]);win<=ceil(wavelengthIn.back());win+=delta)
+ wavelength_fine.push_back(win);
+
+ if(verbose)
+ std::cout << "interpolate wavelengths to " << wavelength_fine.size() << " entries " << std::endl;
+ std::vector<double> srf_fine;//spectral response function, interpolated for wavelength_fine
+
+ stat.interpolateUp(srf[0],srf[1],wavelength_fine,interpolationType,srf_fine,verbose);
+ assert(srf_fine.size()==wavelength_fine.size());
+
+ gsl_interp_accel *accOut;
+ stat.allocAcc(accOut);
+ gsl_spline *splineOut;
+ stat.getSpline(interpolationType,wavelength_fine.size(),splineOut);
+ assert(splineOut);
+
+ std::vector<double> wavelengthOut(wavelength_fine.size());
+
+ for(int iband=0;iband<wavelengthOut.size();++iband)
+ wavelengthOut[iband]=wavelength_fine[iband]*srf_fine[iband];
+
+ stat.initSpline(splineOut,&(wavelength_fine[0]),&(wavelengthOut[0]),wavelength_fine.size());
+ double centreWavelength=gsl_spline_eval_integ(splineOut,start,end,accOut)/norm;
+
+ gsl_spline_free(splineOut);
+ gsl_interp_accel_free(accOut);
+
+ return(centreWavelength);
+}
+
// void filter::Filter::applyFwhm(const vector<double> &wavelengthIn, const ImgReaderGdal& input, const vector<double> &wavelengthOut, const vector<double> &fwhm, const std::string& interpolationType, ImgWriterGdal& output, bool verbose){
// Vector2d<double> lineInput(input.nrOfBand(),input.nrOfCol());
// Vector2d<double> lineOutput(wavelengthOut.size(),input.nrOfCol());
diff --git a/src/algorithms/Filter.h b/src/algorithms/Filter.h
index 26e530f..73bf1ce 100644
--- a/src/algorithms/Filter.h
+++ b/src/algorithms/Filter.h
@@ -61,7 +61,7 @@ public:
template<class T> void morphology(const std::vector<T>& input, std::vector<T>& output, const std::string& method, int dim, short down=1, int offset=0, bool verbose=0);
void morphology(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, int dim, short down=1, int offset=0);
void doit(const ImgReaderGdal& input, ImgWriterGdal& output, short down=1, int offset=0);
-
+ 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);
diff --git a/src/algorithms/StatFactory.h b/src/algorithms/StatFactory.h
index 504c6c7..ad8eef1 100644
--- a/src/algorithms/StatFactory.h
+++ b/src/algorithms/StatFactory.h
@@ -32,9 +32,9 @@ along with pktools. If not, see <http://www.gnu.org/licenses/>.
#include <gsl/gsl_fit.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
-
-using namespace std;
-// using namespace NEWMAT;
+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_randist.h>
+#include <gsl/gsl_cdf.h>
namespace statfactory
{
@@ -43,15 +43,21 @@ class StatFactory{
public:
enum INTERPOLATION_TYPE {undefined=0,linear=1,polynomial=2,cspline=3,cspline_periodic=4,akima=5,akima_periodic=6};
+ //todo: expand with other distributions as in http://www.gnu.org/software/gsl/manual/gsl-ref_toc.html#TOC320
+ enum DISTRIBUTION_TYPE {uniform=1,gaussian=2};
+
StatFactory(void){};
virtual ~StatFactory(void){};
INTERPOLATION_TYPE getInterpolationType(const std::string interpolationType){
std::map<std::string, INTERPOLATION_TYPE> m_interpMap;
initMap(m_interpMap);
- // gsl_interp_accel *acc=gsl_interp_accel_alloc();
- // gsl_spline *spline=gsl_spline_alloc(m_interpMap["interpolationType"],theSize);
return m_interpMap[interpolationType];
};
+ DISTRIBUTION_TYPE getDistributionType(const std::string distributionType){
+ std::map<std::string, DISTRIBUTION_TYPE> m_distMap;
+ initDist(m_distMap);
+ return m_distMap[distributionType];
+ };
static void allocAcc(gsl_interp_accel *&acc){
acc = gsl_interp_accel_alloc ();
};
@@ -89,46 +95,72 @@ public:
return gsl_spline_eval (spline, x, acc);
};
+ static gsl_rng* getRandomGenerator(unsigned long int theSeed){
+ const gsl_rng_type * T;
+ gsl_rng * r;
+
+ // select random number generator
+ r = gsl_rng_alloc (gsl_rng_mt19937);
+ gsl_rng_set(r,theSeed);
+ return r;
+ }
- template<class T> T max(const vector<T>& v) const;
- template<class T> T min(const vector<T>& v) const;
-// template<class T> typename vector<T>::const_iterator max(const vector<T>& v, typename vector<T>::const_iterator begin, typename vector<T>::const_iterator end) const;
- template<class T> typename vector<T>::const_iterator max(const vector<T>& v, typename vector<T>::const_iterator begin, typename vector<T>::const_iterator end) const;
- template<class T> typename vector<T>::iterator max(const vector<T>& v, typename vector<T>::iterator begin, typename vector<T>::iterator end) const;
- template<class T> typename vector<T>::const_iterator absmax(const vector<T>& v, typename vector<T>::const_iterator begin, typename vector<T>::const_iterator end) const;
- template<class T> typename vector<T>::const_iterator min(const vector<T>& v, typename vector<T>::const_iterator begin, typename vector<T>::const_iterator end) const;
- template<class T> typename vector<T>::iterator min(const vector<T>& v, typename vector<T>::iterator begin, typename vector<T>::iterator end) const;
- template<class T> typename vector<T>::const_iterator absmin(const vector<T>& v, typename vector<T>::const_iterator begin, typename vector<T>::const_iterator end) const;
- template<class T> void minmax(const vector<T>& v, typename vector<T>::const_iterator begin, typename vector<T>::const_iterator end, T& theMin, T& theMax) const;
- template<class T> T sum(const vector<T>& v) const;
- template<class T> double mean(const vector<T>& v) const;
- template<class T> T median(const vector<T>& v) const;
- template<class T> double var(const vector<T>& v) const;
- template<class T> double moment(const vector<T>& v, int n) const;
- template<class T> double cmoment(const vector<T>& v, int n) const;
- template<class T> double skewness(const vector<T>& v) const;
- template<class T> double kurtosis(const vector<T>& v) const;
- template<class T> void meanVar(const vector<T>& v, double& m1, double& v1) const;
- template<class T1, class T2> void scale2byte(const vector<T1>& input, vector<T2>& output, unsigned char lbound=0, unsigned char ubound=255) const;
- template<class T> void distribution(const vector<T>& input, typename vector<T>::const_iterator begin, typename vector<T>::const_iterator end, vector<int>& output, int nbin, T &minimum=0.0, T &maximum=0.0, const string &filename="");
- template<class T> void cumulative (const vector<T>& input, typename vector<T>::const_iterator begin, typename vector<T>::const_iterator end, vector<int>& output, int nbin, T &minimum, T &maximum);
- template<class T> void percentiles (const vector<T>& input, typename vector<T>::const_iterator begin, typename vector<T>::const_iterator end, vector<T>& output, int nbin=10, T &minimum=0.0, T &maximum=0.0, const string &filename="");
- template<class T> void signature(const vector<T>& input, double& k, double& alpha, double& beta, double e);
+ static double getRandomValue(const gsl_rng* r, const std::string type, double a=0, double b=0){
+ std::map<std::string, DISTRIBUTION_TYPE> m_distMap;
+ initDist(m_distMap);
+ double randValue=0;
+ switch(m_distMap[type]){
+ case(uniform):
+ randValue = gsl_rng_uniform(r);
+ break;
+ case(gaussian):
+ default:
+ randValue = gsl_ran_gaussian(r, a);
+ break;
+ }
+ return randValue;
+ };
+ template<class T> T max(const std::vector<T>& v) const;
+ template<class T> T min(const std::vector<T>& v) const;
+// template<class T> typename std::vector<T>::const_iterator max(const std::vector<T>& v, typename std::vector<T>::const_iterator begin, typename std::vector<T>::const_iterator end) const;
+ template<class T> typename std::vector<T>::const_iterator max(const std::vector<T>& v, typename std::vector<T>::const_iterator begin, typename std::vector<T>::const_iterator end) const;
+ template<class T> typename std::vector<T>::iterator max(const std::vector<T>& v, typename std::vector<T>::iterator begin, typename std::vector<T>::iterator end) const;
+ template<class T> typename std::vector<T>::const_iterator absmax(const std::vector<T>& v, typename std::vector<T>::const_iterator begin, typename std::vector<T>::const_iterator end) const;
+ template<class T> typename std::vector<T>::const_iterator min(const std::vector<T>& v, typename std::vector<T>::const_iterator begin, typename std::vector<T>::const_iterator end) const;
+ template<class T> typename std::vector<T>::iterator min(const std::vector<T>& v, typename std::vector<T>::iterator begin, typename std::vector<T>::iterator end) const;
+ template<class T> typename std::vector<T>::const_iterator absmin(const std::vector<T>& v, typename std::vector<T>::const_iterator begin, typename std::vector<T>::const_iterator end) const;
+ template<class T> void minmax(const std::vector<T>& v, typename std::vector<T>::const_iterator begin, typename std::vector<T>::const_iterator end, T& theMin, T& theMax) const;
+ template<class T> T sum(const std::vector<T>& v) const;
+ template<class T> double mean(const std::vector<T>& v) const;
+ template<class T> T median(const std::vector<T>& v) const;
+ template<class T> double var(const std::vector<T>& v) const;
+ template<class T> double moment(const std::vector<T>& v, int n) const;
+ template<class T> double cmoment(const std::vector<T>& v, int n) const;
+ template<class T> double skewness(const std::vector<T>& v) const;
+ template<class T> double kurtosis(const std::vector<T>& v) const;
+ template<class T> void meanVar(const std::vector<T>& v, double& m1, double& v1) const;
+ template<class T1, class T2> void scale2byte(const std::vector<T1>& input, std::vector<T2>& output, unsigned char lbound=0, unsigned char ubound=255) const;
+ template<class T> void distribution(const std::vector<T>& input, typename std::vector<T>::const_iterator begin, typename std::vector<T>::const_iterator end, std::vector<double>& output, int nbin, T &minimum=0.0, T &maximum=0.0, double sigma=0, const std::string &filename="");
+ template<class T> void distribution(const std::vector<T>& input, std::vector<double>& output, int nbin, double sigma=0, const std::string &filename=""){distribution(input,input.begin(),input.end(),output,nbin,0,0,sigma,filename);};
+ template<class T> void distribution2d(const std::vector<T>& inputX, const std::vector<T>& inputY, std::vector< std::vector<double> >& output, int nbin, T& minX=0, T& maxX=0, T& minY=0, T& maxY=0, double sigma=0, const std::string& filename="");
+ template<class T> void cumulative (const std::vector<T>& input, typename std::vector<T>::const_iterator begin, typename std::vector<T>::const_iterator end, std::vector<int>& output, int nbin, T &minimum, T &maximum);
+ template<class T> void percentiles (const std::vector<T>& input, typename std::vector<T>::const_iterator begin, typename std::vector<T>::const_iterator end, std::vector<T>& output, int nbin=10, T &minimum=0.0, T &maximum=0.0, const std::string &filename="");
+ template<class T> void signature(const std::vector<T>& input, double& k, double& alpha, double& beta, double e);
void signature(double m1, double m2, double& k, double& alpha, double& beta, double e);
- template<class T> void normalize(const vector<T>& input, vector<double>& output);
- template<class T> void normalize_pct(vector<T>& input);
- template<class T> double rmse(const vector<T>& x, const vector<T>& y) const;
- template<class T> double correlation(const vector<T>& x, const vector<T>& y, int delay=0) const;
- template<class T> double cross_correlation(const vector<T>& x, const vector<T>& y, int maxdelay, vector<T>& z) const;
- template<class T> double linear_regression(const vector<T>& x, const vector<T>& y, double &c0, double &c1) const;
- template<class T> void interpolateUp(const vector<double>& wavelengthIn, const vector<T>& input, const vector<double>& wavelengthOut, const std::string& type, vector<T>& output, bool verbose=false);
- template<class T> void interpolateUp(const vector<double>& wavelengthIn, const vector< vector<T> >& input, const vector<double>& wavelengthOut, const std::string& type, vector< vector<T> >& output, bool verbose=false);
- // template<class T> void interpolateUp(const vector< vector<T> >& input, vector< vector<T> >& output, double start, double end, double step, const gsl_interp_type* type);
- // template<class T> void interpolateUp(const vector< vector<T> >& input, const vector<double>& wavelengthIn, vector< vector<T> >& output, vector<double>& wavelengthOut, double start, double end, double step, const gsl_interp_type* type);
- template<class T> void interpolateUp(const vector<T>& input, vector<T>& output, int nbin);
- template<class T> void interpolateUp(double* input, int dim, vector<T>& output, int nbin);
- template<class T> void interpolateDown(const vector<T>& input, vector<T>& output, int nbin);
- template<class T> void interpolateDown(double* input, int dim, vector<T>& output, int nbin);
+ template<class T> void normalize(const std::vector<T>& input, std::vector<double>& output);
+ template<class T> void normalize_pct(std::vector<T>& input);
+ template<class T> double rmse(const std::vector<T>& x, const std::vector<T>& y) const;
+ template<class T> double correlation(const std::vector<T>& x, const std::vector<T>& y, int delay=0) const;
+ template<class T> double 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> 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);
+ 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);
+ // 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);
+ // template<class T> void interpolateUp(const std::vector< std::vector<T> >& input, const std::vector<double>& wavelengthIn, std::vector< std::vector<T> >& output, std::vector<double>& wavelengthOut, double start, double end, double step, const gsl_interp_type* type);
+ template<class T> void interpolateUp(const std::vector<T>& input, std::vector<T>& output, int nbin);
+ template<class T> void interpolateUp(double* input, int dim, std::vector<T>& output, int nbin);
+ template<class T> void interpolateDown(const std::vector<T>& input, std::vector<T>& output, int nbin);
+ template<class T> void interpolateDown(double* input, int dim, std::vector<T>& output, int nbin);
private:
static void initMap(std::map<std::string, INTERPOLATION_TYPE>& m_interpMap){
@@ -140,74 +172,79 @@ private:
m_interpMap["akima"]=akima;
m_interpMap["akima_periodic"]=akima_periodic;
}
+ static void initDist(std::map<std::string, DISTRIBUTION_TYPE>& m_distMap){
+ //initialize distMap
+ m_distMap["gaussian"]=gaussian;
+ m_distMap["uniform"]=uniform;
+ }
};
-template<class T> typename vector<T>::iterator StatFactory::max(const vector<T>& v, typename vector<T>::iterator begin, typename vector<T>::iterator end) const
+template<class T> typename std::vector<T>::iterator StatFactory::max(const std::vector<T>& v, typename std::vector<T>::iterator begin, typename std::vector<T>::iterator end) const
{
- typename vector<T>::iterator tmpIt=begin;
- for (typename vector<T>::iterator it = begin; it!=end; ++it){
+ typename std::vector<T>::iterator tmpIt=begin;
+ for (typename std::vector<T>::iterator it = begin; it!=end; ++it){
if(*tmpIt<*it)
tmpIt=it;
}
return tmpIt;
}
-template<class T> T StatFactory::max(const vector<T>& v) const
+template<class T> T StatFactory::max(const std::vector<T>& v) const
{
T maxValue=*(v.begin());
- for (typename vector<T>::const_iterator it = v.begin(); it!=v.end(); ++it){
+ for (typename std::vector<T>::const_iterator it = v.begin(); it!=v.end(); ++it){
if(maxValue<*it)
maxValue=*it;
}
return maxValue;
}
-template<class T> T StatFactory::min(const vector<T>& v) const
+template<class T> T StatFactory::min(const std::vector<T>& v) const
{
T minValue=*(v.begin());
- for (typename vector<T>::const_iterator it = v.begin(); it!=v.end(); ++it){
+ for (typename std::vector<T>::const_iterator it = v.begin(); it!=v.end(); ++it){
if(minValue>*it)
minValue=*it;
}
return minValue;
}
-template<class T> typename vector<T>::const_iterator StatFactory::absmax(const vector<T>& v, typename vector<T>::const_iterator begin, typename vector<T>::const_iterator end) const
+template<class T> typename std::vector<T>::const_iterator StatFactory::absmax(const std::vector<T>& v, typename std::vector<T>::const_iterator begin, typename std::vector<T>::const_iterator end) const
{
- typename vector<T>::const_iterator tmpIt=begin;
- for (typename vector<T>::const_iterator it = begin; it!=end; ++it){
+ typename std::vector<T>::const_iterator tmpIt=begin;
+ for (typename std::vector<T>::const_iterator it = begin; it!=end; ++it){
if(abs(*tmpIt)<abs(*it))
tmpIt=it;
}
return tmpIt;
}
-template<class T> typename vector<T>::const_iterator StatFactory::min(const vector<T>& v, typename vector<T>::const_iterator begin, typename vector<T>::const_iterator end) const
+template<class T> typename std::vector<T>::const_iterator StatFactory::min(const std::vector<T>& v, typename std::vector<T>::const_iterator begin, typename std::vector<T>::const_iterator end) const
{
- typename vector<T>::const_iterator tmpIt=begin;
- for (typename vector<T>::const_iterator it = begin; it!=end; ++it){
+ typename std::vector<T>::const_iterator tmpIt=begin;
+ for (typename std::vector<T>::const_iterator it = begin; it!=end; ++it){
if(*tmpIt>*it)
tmpIt=it;
}
return tmpIt;
}
-template<class T> typename vector<T>::const_iterator StatFactory::absmin(const vector<T>& v, typename vector<T>::const_iterator begin, typename vector<T>::const_iterator end) const
+template<class T> typename std::vector<T>::const_iterator StatFactory::absmin(const std::vector<T>& v, typename std::vector<T>::const_iterator begin, typename std::vector<T>::const_iterator end) const
{
- typename vector<T>::const_iterator tmpIt=begin;
- for (typename vector<T>::const_iterator it = begin; it!=end; ++it){
+ typename std::vector<T>::const_iterator tmpIt=begin;
+ for (typename std::vector<T>::const_iterator it = begin; it!=end; ++it){
if(abs(*tmpIt)>abs(*it))
tmpIt=it;
}
}
-template<class T> void StatFactory::minmax(const vector<T>& v, typename vector<T>::const_iterator begin, typename vector<T>::const_iterator end, T& theMin, T& theMax) const
+template<class T> void StatFactory::minmax(const std::vector<T>& v, typename std::vector<T>::const_iterator begin, typename std::vector<T>::const_iterator end, T& theMin, T& theMax) const
{
theMin=*begin;
theMax=*begin;
- for (typename vector<T>::const_iterator it = begin; it!=end; ++it){
+ for (typename std::vector<T>::const_iterator it = begin; it!=end; ++it){
if(theMin>*it)
theMin=*it;
if(theMax<*it)
@@ -215,24 +252,24 @@ template<class T> void StatFactory::minmax(const vector<T>& v, typename vector<T
}
}
-template<class T> T StatFactory::sum(const vector<T>& v) const
+template<class T> T StatFactory::sum(const std::vector<T>& v) const
{
- typename vector<T>::const_iterator it;
+ typename std::vector<T>::const_iterator it;
T tmpSum=0;
for (it = v.begin(); it!= v.end(); ++it)
tmpSum+=*it;
return tmpSum;
}
-template<class T> double StatFactory::mean(const vector<T>& v) const
+template<class T> double StatFactory::mean(const std::vector<T>& v) const
{
assert(v.size());
return static_cast<double>(sum(v))/v.size();
}
-template<class T> T StatFactory::median(const vector<T>& v) const
+template<class T> T StatFactory::median(const std::vector<T>& v) const
{
- vector<T> tmpV=v;
+ std::vector<T> tmpV=v;
sort(tmpV.begin(),tmpV.end());
if(tmpV.size()%2)
return tmpV[tmpV.size()/2];
@@ -240,9 +277,9 @@ template<class T> T StatFactory::median(const vector<T>& v) const
return 0.5*(tmpV[tmpV.size()/2-1]+tmpV[tmpV.size()/2]);
}
-template<class T> double StatFactory::var(const vector<T>& v) const
+template<class T> double StatFactory::var(const std::vector<T>& v) const
{
- typename vector<T>::const_iterator it;
+ typename std::vector<T>::const_iterator it;
double v1=0;
double m1=mean(v);
double n=v.size();
@@ -252,17 +289,17 @@ template<class T> double StatFactory::var(const vector<T>& v) const
v1/=(n-1);
// if(v1<0){
// for (it = v.begin(); it!= v.end(); ++it)
-// cout << *it << " ";
-// cout << endl;
+// std::cout << *it << " ";
+// std::cout << std::endl;
// }
assert(v1>=0);
return v1;
}
-template<class T> double StatFactory::moment(const vector<T>& v, int n) const
+template<class T> double StatFactory::moment(const std::vector<T>& v, int n) const
{
assert(v.size());
- typename vector<T>::const_iterator it;
+ typename std::vector<T>::const_iterator it;
double m=0;
// double m1=mean(v);
for(it = v.begin(); it!= v.end(); ++it){
@@ -272,10 +309,10 @@ template<class T> double StatFactory::moment(const vector<T>& v, int n) const
}
//central moment
-template<class T> double StatFactory::cmoment(const vector<T>& v, int n) const
+template<class T> double StatFactory::cmoment(const std::vector<T>& v, int n) const
{
assert(v.size());
- typename vector<T>::const_iterator it;
+ typename std::vector<T>::const_iterator it;
double m=0;
double m1=mean(v);
for(it = v.begin(); it!= v.end(); ++it){
@@ -284,21 +321,21 @@ template<class T> double StatFactory::cmoment(const vector<T>& v, int n) const
return m/v.size();
}
-template<class T> double StatFactory::skewness(const vector<T>& v) const
+template<class T> double StatFactory::skewness(const std::vector<T>& v) const
{
return cmoment(v,3)/pow(var(v),1.5);
}
-template<class T> double StatFactory::kurtosis(const vector<T>& v) const
+template<class T> double StatFactory::kurtosis(const std::vector<T>& v) const
{
double m2=cmoment(v,2);
double m4=cmoment(v,4);
return m4/m2/m2-3.0;
}
-template<class T> void StatFactory::meanVar(const vector<T>& v, double& m1, double& v1) const
+template<class T> void StatFactory::meanVar(const std::vector<T>& v, double& m1, double& v1) const
{
- typename vector<T>::const_iterator it;
+ typename std::vector<T>::const_iterator it;
v1=0;
m1=mean(v);
double n=v.size();
@@ -309,7 +346,7 @@ template<class T> void StatFactory::meanVar(const vector<T>& v, double& m1, doub
assert(v1>=0);
}
-template<class T1, class T2> void StatFactory::scale2byte(const vector<T1>& input, vector<T2>& output, unsigned char lbound, unsigned char ubound) const
+template<class T1, class T2> void StatFactory::scale2byte(const std::vector<T1>& input, std::vector<T2>& output, unsigned char lbound, unsigned char ubound) const
{
output.resize(input.size());
T1 minimum=min(input);
@@ -320,7 +357,7 @@ template<class T1, class T2> void StatFactory::scale2byte(const vector<T1>& inpu
output[i]=scale*(input[i]-(minimum))+lbound;
}
-template<class T> void StatFactory::distribution (const vector<T>& input, typename vector<T>::const_iterator begin, typename vector<T>::const_iterator end, vector<int>& output, int nbin, T &minimum, T &maximum, const string &filename)
+template<class T> void StatFactory::distribution(const std::vector<T>& input, typename std::vector<T>::const_iterator begin, typename std::vector<T>::const_iterator end, std::vector<double>& output, int nbin, T &minimum, T &maximum, double sigma, const std::string &filename)
{
if(maximum<=minimum)
minmax(input,begin,end,minimum,maximum);
@@ -329,8 +366,8 @@ template<class T> void StatFactory::distribution (const vector<T>& input, typen
// if(!maximum)
// maximum=*(max(input,begin,end));
if(maximum<=minimum){
- ostringstream s;
- s<<"Error: could not calculate distribution (min>=max) in " << filename;
+ std::ostringstream s;
+ s<<"Error: could not calculate distribution (min>=max)";
throw(s.str());
}
assert(nbin>1);
@@ -339,28 +376,128 @@ template<class T> void StatFactory::distribution (const vector<T>& input, typen
output.resize(nbin);
for(int i=0;i<nbin;output[i++]=0);
}
- typename vector<T>::const_iterator it;
+ typename std::vector<T>::const_iterator it;
for(it=begin;it!=end;++it){
- if(*it==maximum)
- ++output[nbin-1];
- else if(*it>=minimum && *it<maximum)
- ++output[static_cast<int>(static_cast<double>((*it)-minimum)/(maximum-minimum)*nbin)];
+ if(sigma>0){
+ // minimum-=2*sigma;
+ // maximum+=2*sigma;
+ //create kde for Gaussian basis function
+ //todo: speed up by calculating first and last bin with non-zero contriubtion...
+ //todo: calculate real surface below pdf by using gsl_cdf_gaussian_P(x-mean+binsize,sigma)-gsl_cdf_gaussian_P(x-mean,sigma)
+ for(int ibin=0;ibin<nbin;++ibin){
+ double icenter=minimum+static_cast<double>(maximum-minimum)*(ibin+0.5)/nbin;
+ double thePdf=gsl_ran_gaussian_pdf(*it-icenter, sigma);
+ output[ibin]+=thePdf;
+ }
+ }
+ else{
+ int theBin=0;
+ if(*it==maximum)
+ theBin=nbin-1;
+ else if(*it>=minimum && *it<maximum)
+ theBin=static_cast<int>(static_cast<double>((*it)-minimum)/(maximum-minimum)*nbin);
+ ++output[theBin];
+ // if(*it==maximum)
+ // ++output[nbin-1];
+ // else if(*it>=minimum && *it<maximum)
+ // ++output[static_cast<int>(static_cast<double>((*it)-minimum)/(maximum-minimum)*nbin)];
+ }
}
if(!filename.empty()){
- ofstream outputfile;
+ std::ofstream outputfile;
outputfile.open(filename.c_str());
if(!outputfile){
- ostringstream s;
+ std::ostringstream s;
s<<"Error opening distribution file , " << filename;
throw(s.str());
}
for(int bin=0;bin<nbin;++bin)
- outputfile << (maximum-minimum)*bin/(nbin-1)+minimum << " " << static_cast<double>(output[bin])/input.size() << endl;
+ outputfile << (maximum-minimum)*bin/(nbin-1)+minimum << " " << static_cast<double>(output[bin])/input.size() << std::endl;
+ outputfile.close();
+ }
+}
+
+template<class T> void StatFactory::distribution2d(const std::vector<T>& inputX, const std::vector<T>& inputY, std::vector< std::vector<double> >& output, int nbin, T& minX, T& maxX, T& minY, T& maxY, double sigma, const std::string& filename)
+{
+ assert(inputX.size());
+ assert(inputX.size()==inputY.size());
+ unsigned long int npoint=inputX.size();
+ if(maxX<=minX)
+ minmax(inputX,inputX.begin(),inputX.end(),minX,maxX);
+ if(maxX<=minX){
+ std::ostringstream s;
+ s<<"Error: could not calculate distribution (minX>=maxX)";
+ throw(s.str());
+ }
+ if(maxY<=minY)
+ minmax(inputY,inputY.begin(),inputY.end(),minY,maxY);
+ if(maxY<=minY){
+ std::ostringstream s;
+ s<<"Error: could not calculate distribution (minY>=maxY)";
+ throw(s.str());
+ }
+ assert(nbin>1);
+ output.resize(nbin);
+ for(int i=0;i<nbin;++i){
+ output[i].resize(nbin);
+ for(int j=0;j<nbin;++j)
+ output[i][j]=0;
+ }
+ int binX=0;
+ int binY=0;
+ for(unsigned long int ipoint=0;ipoint<npoint;++ipoint){
+ if(inputX[ipoint]==maxX)
+ binX=nbin-1;
+ else
+ binX=static_cast<int>(static_cast<double>(inputX[ipoint]-minX)/(maxX-minX)*nbin);
+ if(inputY[ipoint]==maxY)
+ binY=nbin-1;
+ else
+ binY=static_cast<int>(static_cast<double>(inputY[ipoint]-minY)/(maxY-minY)*nbin);
+ assert(binX>=0);
+ assert(binX<output.size());
+ assert(binY>=0);
+ assert(binY<output[binX].size());
+ if(sigma>0){
+ // minX-=2*sigma;
+ // maxX+=2*sigma;
+ // minY-=2*sigma;
+ // maxY+=2*sigma;
+ //create kde for Gaussian basis function
+ //todo: speed up by calculating first and last bin with non-zero contriubtion...
+ for(int ibinX=0;ibinX<nbin;++ibinX){
+ double centerX=minX+static_cast<double>(maxX-minX)*ibinX/nbin;
+ double pdfX=gsl_ran_gaussian_pdf(inputX[ipoint]-centerX, sigma);
+ for(int ibinY=0;ibinY<nbin;++ibinY){
+ //calculate \integral_ibinX^(ibinX+1)
+ double centerY=minY+static_cast<double>(maxY-minY)*ibinY/nbin;
+ double pdfY=gsl_ran_gaussian_pdf(inputY[ipoint]-centerY, sigma);
+ output[ibinX][binY]+=pdfX*pdfY;
+ }
+ }
+ }
+ else
+ ++output[binX][binY];
+ }
+ if(!filename.empty()){
+ std::ofstream outputfile;
+ outputfile.open(filename.c_str());
+ if(!outputfile){
+ std::ostringstream s;
+ s<<"Error opening distribution file , " << filename;
+ throw(s.str());
+ }
+ for(int bin=0;bin<nbin;++bin){
+ for(int bin=0;bin<nbin;++bin){
+ double value=static_cast<double>(output[binX][binY])/npoint;
+ outputfile << (maxX-minX)*bin/(nbin-1)+minX << " " << (maxY-minY)*bin/(nbin-1)+minY << " " << value << std::endl;
+ }
+ }
outputfile.close();
}
}
-template<class T> void StatFactory::percentiles (const vector<T>& input, typename vector<T>::const_iterator begin, typename vector<T>::const_iterator end, vector<T>& output, int nbin, T &minimum, T &maximum, const string &filename)
+template<class T> void StatFactory::percentiles (const std::vector<T>& input, typename std::vector<T>::const_iterator begin, typename std::vector<T>::const_iterator end, std::vector<T>& output, int nbin, T &minimum, T &maximum, const std::string &filename)
{
if(maximum<=minimum)
minmax(input,begin,end,minimum,maximum);
@@ -394,48 +531,48 @@ template<class T> void StatFactory::percentiles (const vector<T>& input, typena
output[ibin]=median(inputBin);
}
if(!filename.empty()){
- ofstream outputfile;
+ std::ofstream outputfile;
outputfile.open(filename.c_str());
if(!outputfile){
- ostringstream s;
+ std::ostringstream s;
s<<"error opening distribution file , " << filename;
throw(s.str());
}
for(int ibin=0;ibin<nbin;++ibin)
- outputfile << ibin*100.0/nbin << " " << static_cast<double>(output[ibin])/input.size() << endl;
+ outputfile << ibin*100.0/nbin << " " << static_cast<double>(output[ibin])/input.size() << std::endl;
outputfile.close();
}
}
-// template<class T> void StatFactory::cumulative (const vector<T>& input, typename vector<T>::const_iterator begin, typename vector<T>::const_iterator end, vector<int>& output, int nbin, T &minimum, T &maximum)
+// template<class T> void StatFactory::cumulative (const std::vector<T>& input, typename std::vector<T>::const_iterator begin, typename std::vector<T>::const_iterator end, std::vector<int>& output, int nbin, T &minimum, T &maximum)
// {
// assert(nbin>1);
// assert(input.size());
// distribution(input,output,nbin,minimum,maximum);
-// for(vector<int>::iterator it=begin+1;it!=end;++it)
+// for(std::vector<int>::iterator it=begin+1;it!=end;++it)
// *it+=*(it-1);
// if(!filename.empty()){
-// ofstream outputfile;
+// std::ofstream outputfile;
// outputfile.open(filename.c_str());
// if(!outputfile){
-// ostringstream s;
+// std::ostringstream s;
// s<<"error opening cumulative file , " << filename;
// throw(s.str());
// }
// for(int bin=0;bin<nbin;++bin)
-// outputfile << (maximum-minimum)*bin/(nbin-1)+minimum << " " << static_cast<double>(output[bin])/input.size() << endl;
+// outputfile << (maximum-minimum)*bin/(nbin-1)+minimum << " " << static_cast<double>(output[bin])/input.size() << std::endl;
// outputfile.close();
// }
// }
-template<class T> void StatFactory::signature(const vector<T>& input, double&k, double& alpha, double& beta, double e)
+template<class T> void StatFactory::signature(const std::vector<T>& input, double&k, double& alpha, double& beta, double e)
{
double m1=moment(input,1);
double m2=moment(input,2);
signature(m1,m2,k,alpha,beta,e);
}
-template<class T> void StatFactory::normalize(const vector<T>& input, vector<double>& output){
+template<class T> void StatFactory::normalize(const std::vector<T>& input, std::vector<double>& output){
double total=sum(input);
if(total){
output.resize(input.size());
@@ -446,16 +583,16 @@ template<class T> void StatFactory::normalize(const vector<T>& input, vector<dou
output=input;
}
-template<class T> void StatFactory::normalize_pct(vector<T>& input){
+template<class T> void StatFactory::normalize_pct(std::vector<T>& input){
double total=sum(input);
if(total){
- typename vector<T>::iterator it;
+ typename std::vector<T>::iterator it;
for(it=input.begin();it!=input.end();++it)
*it=100.0*(*it)/total;
}
}
-template<class T> double StatFactory::rmse(const vector<T>& x, const vector<T>& y) const{
+template<class T> double StatFactory::rmse(const std::vector<T>& x, const std::vector<T>& y) const{
assert(x.size()==y.size());
assert(x.size());
double mse=0;
@@ -466,7 +603,7 @@ template<class T> double StatFactory::rmse(const vector<T>& x, const vector<T>&
return sqrt(mse);
}
-template<class T> double StatFactory::correlation(const vector<T>& x, const vector<T>& y, int delay) const{
+template<class T> double StatFactory::correlation(const std::vector<T>& x, const std::vector<T>& y, int delay) const{
double meanX=0;
double meanY=0;
double varX=0;
@@ -495,7 +632,7 @@ template<class T> double StatFactory::correlation(const vector<T>& x, const vect
return 0;
}
-template<class T> double StatFactory::cross_correlation(const vector<T>& x, const vector<T>& y, int maxdelay, vector<T>& z) const{
+template<class T> double StatFactory::cross_correlation(const std::vector<T>& x, const std::vector<T>& y, int maxdelay, std::vector<T>& z) const{
z.clear();
double sumCorrelation=0;
for (int delay=-maxdelay;delay<maxdelay;delay++) {
@@ -505,7 +642,7 @@ template<class T> double StatFactory::cross_correlation(const vector<T>& x, cons
return sumCorrelation;
}
- template<class T> double StatFactory::linear_regression(const vector<T>& x, const 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;
@@ -519,7 +656,7 @@ template<class T> double StatFactory::cross_correlation(const vector<T>& x, cons
//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)
-template<class T> void StatFactory::interpolateUp(const vector<double>& wavelengthIn, const vector<T>& input, const vector<double>& wavelengthOut, const std::string& type, vector<T>& output, bool verbose){
+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){
assert(wavelengthIn.size());
assert(input.size()==wavelengthIn.size());
assert(wavelengthOut.size());
@@ -551,7 +688,7 @@ template<class T> void StatFactory::interpolateUp(const vector<double>& waveleng
gsl_interp_accel_free(acc);
}
-// template<class T> void StatFactory::interpolateUp(const vector<double>& wavelengthIn, const vector< vector<T> >& input, const vector<double>& wavelengthOut, const std::string& type, vector< vector<T> >& output, bool verbose){
+// template<class T> void StatFactory::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){
// assert(wavelengthIn.size());
// assert(wavelengthOut.size());
// int nsample=input.size();
@@ -582,7 +719,7 @@ template<class T> void StatFactory::interpolateUp(const vector<double>& waveleng
// gsl_interp_accel_free(acc);
// }
-template<class T> void StatFactory::interpolateUp(const vector<T>& input, vector<T>& output, int nbin)
+template<class T> void StatFactory::interpolateUp(const std::vector<T>& input, std::vector<T>& output, int nbin)
{
assert(input.size());
assert(nbin);
@@ -603,7 +740,7 @@ template<class T> void StatFactory::interpolateUp(const vector<T>& input, vector
}
}
-template<class T> void StatFactory::interpolateUp(double* input, int dim, vector<T>& output, int nbin)
+template<class T> void StatFactory::interpolateUp(double* input, int dim, std::vector<T>& output, int nbin)
{
assert(nbin);
output.clear();
@@ -622,7 +759,7 @@ template<class T> void StatFactory::interpolateUp(double* input, int dim, vector
}
}
-template<class T> void StatFactory::interpolateDown(const vector<T>& input, vector<T>& output, int nbin)
+template<class T> void StatFactory::interpolateDown(const std::vector<T>& input, std::vector<T>& output, int nbin)
{
assert(input.size());
assert(nbin);
@@ -640,7 +777,7 @@ template<class T> void StatFactory::interpolateDown(const vector<T>& input, vect
}
}
-template<class T> void StatFactory::interpolateDown(double* input, int dim, vector<T>& output, int nbin)
+template<class T> void StatFactory::interpolateDown(double* input, int dim, std::vector<T>& output, int nbin)
{
assert(nbin);
output.clear();
@@ -667,7 +804,7 @@ template<class T> void StatFactory::interpolateDown(double* input, int dim, vect
// double g=exp(lgamma(1.0/beta));
// alpha=m1*g/exp(lgamma(2.0/beta));
// k=beta/(2*alpha*g);
-// // cout << "y, alpha, beta: " << y << ", " << alpha << ", " << beta << endl;
+// // std::cout << "y, alpha, beta: " << y << ", " << alpha << ", " << beta << std::endl;
// }
// double Histogram::F(double x)
diff --git a/src/algorithms/myfann_cpp.h b/src/algorithms/myfann_cpp.h
index 1308e22..e3d4005 100644
--- a/src/algorithms/myfann_cpp.h
+++ b/src/algorithms/myfann_cpp.h
@@ -1540,7 +1540,7 @@ public:
assert(cv<ntraining);
float rmse=0;
int nclass=trainingFeatures.size();
- vector< Vector2d<float> > testFeatures(nclass);
+ std::vector< Vector2d<float> > testFeatures(nclass);
int testclass=0;//class to leave out
int testsample=0;//sample to leave out
int nrun=(cv>1)? cv : ntraining;
@@ -1589,16 +1589,16 @@ public:
assert(nsample==ntest);
//training with left out training set
if(verbose>1)
- cout << endl << "Set training data" << endl;
+ std::cout << std::endl << "Set training data" << std::endl;
bool initWeights=true;
unsigned int epochs_between_reports=0;
train_on_data(trainingFeatures,ntraining-ntest,initWeights, max_epochs,
epochs_between_reports, desired_error);
//cross validation with testFeatures
if(verbose>1)
- cout << endl << "Cross validation" << endl;
+ std::cout << std::endl << "Cross validation" << std::endl;
- vector<float> result(nclass);
+ std::vector<float> result(nclass);
int maxClass=-1;
for(int iclass=0;iclass<testFeatures.size();++iclass){
assert(trainingFeatures[iclass].size());
@@ -1621,7 +1621,7 @@ public:
// rmse+=test_data(testFeatures,ntest);
// if(verbose>1)
- // cout << endl << "rmse: " << rmse << endl;
+ // std::cout << std::endl << "rmse: " << rmse << std::endl;
}
// rmse/=nrun;
//reset from very last run
@@ -1698,7 +1698,7 @@ public:
assert(testInput.size()==testOutput.size());
//training with left out training set
if(verbose>1)
- cout << endl << "Set training data" << endl;
+ std::cout << std::endl << "Set training data" << std::endl;
bool initWeights=true;
unsigned int epochs_between_reports=0;
@@ -1706,9 +1706,9 @@ public:
epochs_between_reports, desired_error);
//cross validation with testFeatures
if(verbose>1)
- cout << endl << "Cross validation" << endl;
+ std::cout << std::endl << "Cross validation" << std::endl;
- vector<fann_type> result(noutput);
+ std::vector<fann_type> result(noutput);
for(int isample=0;isample<testInput.size();++isample){
result=run(testInput[isample]);
referenceVector.push_back(testOutput[isample]);
diff --git a/src/apps/Makefile.am b/src/apps/Makefile.am
index 91e1db9..eeafcff 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 pkstat pkstatogr pkegcs pkextract pkfillnodata pkfilter pkenhance pkfilterascii pkdsm2shadow pkmosaic pkndvi pkpolygonize pkascii2img pkdiff pkclassify_svm pkfs_svm pkascii2ogr pkeditogr
+bin_PROGRAMS = pkinfo pkcrop pkreclass pkgetmask pksetmask pkcreatect pkdumpimg pkdumpogr pksieve pkstatascii pkstatogr pkegcs pkextract pkfillnodata pkfilter pkenhance pkfilterascii pkdsm2shadow pkmosaic pkndvi pkpolygonize pkascii2img pkdiff pkclassify_svm pkfs_svm pkascii2ogr pkeditogr
# the program to build but not install (the names of the final binaries)
#noinst_PROGRAMS = pkxcorimg pkgeom
@@ -44,8 +44,8 @@ pkcreatect_SOURCES = pkcreatect.cc
pkdumpimg_SOURCES = pkdumpimg.cc
pkdumpogr_SOURCES = pkdumpogr.h pkdumpogr.cc
pksieve_SOURCES = pksieve.cc
-pkstat_SOURCES = $(top_srcdir)/src/algorithms/StatFactory.h pkstat.cc
-pkstat_LDADD = $(GSL_LIBS) $(AM_LDFLAGS)
+pkstatascii_SOURCES = $(top_srcdir)/src/algorithms/StatFactory.h pkstatascii.cc
+pkstatascii_LDADD = $(GSL_LIBS) $(AM_LDFLAGS)
pkstatogr_SOURCES = pkstatogr.cc
pkegcs_SOURCES = pkegcs.cc
pkegcs_LDADD = -lgdal $(AM_LDFLAGS) -lgdal
diff --git a/src/apps/pkegcs.cc b/src/apps/pkegcs.cc
index 79f362b..6129a6c 100644
--- a/src/apps/pkegcs.cc
+++ b/src/apps/pkegcs.cc
@@ -23,10 +23,10 @@ along with pktools. If not, see <http://www.gnu.org/licenses/>.
int main(int argc, char *argv[])
{
- Optionpk<string> image_opt("i","image","input image to analyse","");
+ Optionpk<std::string> image_opt("i","image","input image to analyse","");
Optionpk<unsigned short> band_opt("b", "band", "Band specific information", 0);
- Optionpk<string> cell2bb_opt("c2b","cell2bb","convert cell code to geo coordinates of boundingbox (e.g. 32-AB)","");
- Optionpk<string> cell2mid_opt("c2m","cell2mid","convert cell code to centre in geo coordinates (e.g. 32-AB)","");
+ Optionpk<std::string> cell2bb_opt("c2b","cell2bb","convert cell code to geo coordinates of boundingbox (e.g. 32-AB)","");
+ Optionpk<std::string> cell2mid_opt("c2m","cell2mid","convert cell code to centre in geo coordinates (e.g. 32-AB)","");
Optionpk<bool> refpixel_opt("\0", "ref", "get reference pixel (lower left corner of centre of gravity pixel)", false);
Optionpk<double> maskValue_opt("m", "mask", "mask value(s) for no data to calculate reference pixel in image",0);
Optionpk<int> dx_opt("dx","dx","resolution",250);
@@ -48,7 +48,7 @@ int main(int argc, char *argv[])
x_opt.retrieveOption(argc,argv);
y_opt.retrieveOption(argc,argv);
}
- catch(string predefinedString){
+ catch(std::string predefinedString){
std::cout << predefinedString << std::endl;
exit(0);
}
@@ -62,17 +62,17 @@ int main(int argc, char *argv[])
int theULX, theULY, theLRX, theLRY;
egcs.setLevel(egcs.cell2level(cell2bb_opt[0]));
egcs.cell2bb(cell2bb_opt[0],theULX,theULY,theLRX,theLRY);
- cout << setprecision(12) << "--ulx=" << theULX << " --uly=" << theULY << " --lrx=" << theLRX << " --lry=" << theLRY << endl;
+ std::cout << std::setprecision(12) << "--ulx=" << theULX << " --uly=" << theULY << " --lrx=" << theLRX << " --lry=" << theLRY << std::endl;
}
if(cell2mid_opt[0]!=""){
double midX, midY;
egcs.setLevel(egcs.cell2level(cell2mid_opt[0]));
egcs.cell2mid(cell2mid_opt[0],midX,midY);
- cout << setprecision(12) << "-x=" << midX << " -y=" << midY << endl;
+ std::cout << std::setprecision(12) << "-x=" << midX << " -y=" << midY << std::endl;
}
if(geo2cell_opt[0]){
egcs.setLevel(egcs.res2level(dx_opt[0]));
- cout << egcs.geo2cell(x_opt[0],y_opt[0]) << endl;
+ std::cout << egcs.geo2cell(x_opt[0],y_opt[0]) << std::endl;
}
if(image_opt[0]!=""){
ImgReaderGdal imgReader;
@@ -84,16 +84,16 @@ int main(int argc, char *argv[])
// if(verbose_opt[0]){
// vector<double> noData;
// imgReader.getNoDataValues(noData,band_opt[0]);
- // cout << "number of no data values: " << noData.size() << endl;
+ // std::cout << "number of no data values: " << noData.size() << std::endl;
// }
double refX,refY;
//get centre of reference (centre of gravity) pixel in image
imgReader.getRefPix(refX,refY,band_opt[0]);
- cout << setprecision(12) << "--x " << refX << " --y " << refY << endl;
+ std::cout << std::setprecision(12) << "--x " << refX << " --y " << refY << std::endl;
egcs.setLevel(egcs.res2level(imgReader.getDeltaX()));
// unsigned short theLevel=egcs.getLevel(imgReader.getDeltaX());
// egcs.setLevel(theLevel);
- cout << "cell code at level " << egcs.getLevel() << " (resolution is " << egcs.getResolution() << "): " << egcs.geo2cell(refX,refY) << endl;
+ std::cout << "cell code at level " << egcs.getLevel() << " (resolution is " << egcs.getResolution() << "): " << egcs.geo2cell(refX,refY) << std::endl;
}
}
}
diff --git a/src/apps/pkinfo.cc b/src/apps/pkinfo.cc
index 6095ad8..b53e89f 100644
--- a/src/apps/pkinfo.cc
+++ b/src/apps/pkinfo.cc
@@ -104,7 +104,7 @@ int main(int argc, char *argv[])
metadata_opt.retrieveOption(argc,argv);
nodata_opt.retrieveOption(argc,argv);
}
- catch(string predefinedString){
+ catch(std::string predefinedString){
std::cout << predefinedString << std::endl;
exit(0);
}
@@ -137,7 +137,7 @@ int main(int argc, char *argv[])
if(centre_opt[0]){
double theX, theY;
imgReader.getCentrePos(theX,theY);
- std::cout << setprecision(12) << " -x " << theX << " -y " << theY << " ";
+ std::cout << std::setprecision(12) << " -x " << theX << " -y " << theY << " ";
}
if(refpixel_opt[0]){
assert(band_opt[0]<imgReader.nrOfBand());
@@ -145,7 +145,7 @@ int main(int argc, char *argv[])
double refX,refY;
//get centre of reference (centre of gravity) pixel in image
imgReader.getRefPix(refX,refY,band_opt[0]);
- cout << setprecision(12) << "-rx " << refX << " -ry " << refY << endl;
+ std::cout << std::setprecision(12) << "-rx " << refX << " -ry " << refY << std::endl;
egcs.setLevel(egcs.res2level(imgReader.getDeltaX()));
// unsigned short theLevel=egcs.getLevel(imgReader.getDeltaX());
// egcs.setLevel(theLevel);
@@ -155,9 +155,9 @@ int main(int argc, char *argv[])
double theULX, theULY, theLRX, theLRY;
imgReader.getBoundingBox(theULX,theULY,theLRX,theLRY);
if(bbox_te_opt[0])
- std::cout << setprecision(12) << "-te " << theULX << " " << theLRY << " " << theLRX << " " << theULY;
+ std::cout << std::setprecision(12) << "-te " << theULX << " " << theLRY << " " << theLRX << " " << theULY;
else
- std::cout << setprecision(12) << "--ulx=" << theULX << " --uly=" << theULY << " --lrx=" << theLRX << " --lry=" << theLRY << " ";
+ std::cout << std::setprecision(12) << "--ulx=" << theULX << " --uly=" << theULY << " --lrx=" << theLRX << " --lry=" << theLRY << " ";
if(!ifile){
maxLRX=theLRX;
maxULY=theULY;
@@ -191,7 +191,7 @@ int main(int argc, char *argv[])
if(extent_opt.size()){
extentReader.open(extent_opt[0]);
if(!(extentReader.getExtent(ulx_opt[0],uly_opt[0],lrx_opt[0],lry_opt[0]))){
- cerr << "Error: could not get extent from " << extent_opt[0] << std::endl;
+ std::cerr << "Error: could not get extent from " << extent_opt[0] << std::endl;
exit(1);
}
// std::cout << "--ulx=" << ulx_opt[0] << " --uly=" << uly_opt[0] << " --lrx=" << lrx_opt[0] << " --lry=" << lry_opt[0] << std::endl;
@@ -208,13 +208,13 @@ int main(int argc, char *argv[])
double ulx,uly,lrx,lry;
imgReader.getBoundingBox(ulx,uly,lrx,lry);
if(ulx_opt.size())
- std::cout << " --ulx=" << fixed << ulx << " ";
+ std::cout << " --ulx=" << std::fixed << ulx << " ";
if(uly_opt.size())
- std::cout << " --uly=" << fixed << uly << " ";
+ std::cout << " --uly=" << std::fixed << uly << " ";
if(lrx_opt.size())
- std::cout << " --lrx=" << fixed << lrx << " ";
+ std::cout << " --lrx=" << std::fixed << lrx << " ";
if(lry_opt.size())
- std::cout << " --lry=" << fixed << lry << " ";
+ std::cout << " --lry=" << std::fixed << lry << " ";
}
if(colorTable_opt[0]){
GDALColorTable* colorTable=imgReader.getColorTable();
@@ -334,7 +334,7 @@ int main(int argc, char *argv[])
if(geo_opt[0]&&!read_opt[0]){
double ulx,uly,deltaX,deltaY,rot1,rot2;
imgReader.getGeoTransform(ulx,uly,deltaX,deltaY,rot1,rot2);
- std::cout << " --geo " << setprecision(12) << ulx << " " << uly << " " << deltaX << " " << deltaY << " " << rot1 << " " << rot2 << " ";
+ std::cout << " --geo " << std::setprecision(12) << ulx << " " << uly << " " << deltaX << " " << deltaY << " " << rot1 << " " << rot2 << " ";
}
if(interleave_opt[0]){
std::cout << " --interleave " << imgReader.getInterleave() << " ";
@@ -350,18 +350,18 @@ int main(int argc, char *argv[])
// catch(...){
// std::cout << "catched" << std::endl;
// }
- list<std::string> metaData;
+ std::list<std::string> metaData;
imgReader.getMetadata(metaData);
- list<std::string>::const_iterator lit=metaData.begin();
+ std::list<std::string>::const_iterator lit=metaData.begin();
std::cout << " --description ";
while(lit!=metaData.end())
std::cout << *(lit++) << " ";
}
if(metadata_opt[0]){
std::cout << "Metadata: " << std::endl;
- list<std::string> lmeta;
+ std::list<std::string> lmeta;
imgReader.getMetadata(lmeta);
- list<std::string>::const_iterator lit=lmeta.begin();
+ std::list<std::string>::const_iterator lit=lmeta.begin();
while(lit!=lmeta.end()){
std::cout << *lit << std::endl;
++lit;
@@ -407,14 +407,14 @@ int main(int argc, char *argv[])
}
if((bbox_opt[0]||bbox_te_opt[0])&&input_opt.size()>1){
if(bbox_te_opt[0])
- std::cout << setprecision(12) << "-te " << minULX << " " << minLRY << " " << maxLRX << " " << maxULY;
+ std::cout << std::setprecision(12) << "-te " << minULX << " " << minLRY << " " << maxLRX << " " << maxULY;
else
- std::cout << "union bounding box: " << setprecision(12) << "--ulx=" << minULX << " --uly=" << maxULY << " --lrx=" << maxLRX << " --lry=" << minLRY << std::endl;
+ std::cout << "union bounding box: " << std::setprecision(12) << "--ulx=" << minULX << " --uly=" << maxULY << " --lrx=" << maxLRX << " --lry=" << minLRY << std::endl;
if(maxULX<minLRX&&minULY>maxLRY){
if(bbox_te_opt[0])
- std::cout << "intersect bounding box: " << setprecision(12) << "-te " << maxULX << " " << maxLRY << " " << minLRX << " --lry=" << minULY << std::endl;
+ std::cout << "intersect bounding box: " << std::setprecision(12) << "-te " << maxULX << " " << maxLRY << " " << minLRX << " --lry=" << minULY << std::endl;
else
- std::cout << "intersect bounding box: " << setprecision(12) << "--ulx=" << maxULX << " --uly=" << minULY << " --lrx=" << minLRX << " --lry=" << maxLRY << std::endl;
+ std::cout << "intersect bounding box: " << std::setprecision(12) << "--ulx=" << maxULX << " --uly=" << minULY << " --lrx=" << minLRX << " --lry=" << maxLRY << std::endl;
}
else
std::cout << "no intersect" << std::endl;
diff --git a/src/apps/pkregression_nn.cc b/src/apps/pkregression_nn.cc
index 35e0a53..4f3192c 100644
--- a/src/apps/pkregression_nn.cc
+++ b/src/apps/pkregression_nn.cc
@@ -24,6 +24,8 @@ along with pktools. If not, see <http://www.gnu.org/licenses/>.
#include "floatfann.h"
#include "myfann_cpp.h"
+using namespace std;
+
int main(int argc, char *argv[])
{
//--------------------------- command line options ------------------------------------
diff --git a/src/apps/pkstat.cc b/src/apps/pkstatascii.cc
similarity index 65%
rename from src/apps/pkstat.cc
rename to src/apps/pkstatascii.cc
index 10caf0d..d5cf526 100644
--- a/src/apps/pkstat.cc
+++ b/src/apps/pkstatascii.cc
@@ -24,7 +24,6 @@ along with pktools. If not, see <http://www.gnu.org/licenses/>.
#include "base/Optionpk.h"
#include "fileclasses/FileReaderAscii.h"
#include "algorithms/StatFactory.h"
-
using namespace std;
int main(int argc, char *argv[])
@@ -36,7 +35,11 @@ int main(int argc, char *argv[])
Optionpk<int> col_opt("c", "column", "column nr, starting from 0", 0);
Optionpk<int> range_opt("r", "range", "rows to start/end reading. Use -r 1 -r 10 to read first 10 rows where first row is header. Use 0 to read all rows with no header.", 0);
Optionpk<bool> size_opt("size","size","sample size",false);
- Optionpk<bool> mean_opt("m","mean","calculate mean value",false);
+ Optionpk<unsigned int> rand_opt("rnd", "rnd", "generate random numbers", 0);
+ Optionpk<std::string> randdist_opt("dist", "dist", "distribution for generating random numbers, see http://www.gn/software/gsl/manual/gsl-ref_toc.html#TOC320 (only uniform and Gaussian supported yet)", "gaussian");
+ Optionpk<double> randa_opt("rnda", "rnda", "first parameter for random distribution (standard deviation in case of Gaussian)", 1);
+ Optionpk<double> randb_opt("rndb", "rndb", "second parameter for random distribution (standard deviation in case of Gaussian)", 0);
+ Optionpk<bool> mean_opt("m","mean","calculate median",false);
Optionpk<bool> median_opt("med","median","calculate median",false);
Optionpk<bool> var_opt("var","var","calculate variance",false);
Optionpk<bool> skewness_opt("skew","skewness","calculate skewness",false);
@@ -47,8 +50,10 @@ int main(int argc, char *argv[])
Optionpk<double> min_opt("min","min","set minimum value",0);
Optionpk<double> max_opt("max","max","set maximum value",0);
Optionpk<bool> histogram_opt("hist","hist","calculate histogram",false);
+ Optionpk<bool> histogram2d_opt("hist2d","hist2d","calculate 2-dimensional histogram based on two columns",false);
Optionpk<short> nbin_opt("bin","bin","number of bins to calculate histogram",0);
Optionpk<bool> relative_opt("rel","relative","use percentiles for histogram to calculate histogram",false);
+ Optionpk<double> kde_opt("kde","kde","bandwith of kernel density when producing histogram, use 0 for practical estimation based on Silverman's rule of thumb");
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("e","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);
@@ -63,6 +68,10 @@ int main(int argc, char *argv[])
col_opt.retrieveOption(argc,argv);
range_opt.retrieveOption(argc,argv);
size_opt.retrieveOption(argc,argv);
+ rand_opt.retrieveOption(argc,argv);
+ randdist_opt.retrieveOption(argc,argv);
+ randa_opt.retrieveOption(argc,argv);
+ randb_opt.retrieveOption(argc,argv);
mean_opt.retrieveOption(argc,argv);
median_opt.retrieveOption(argc,argv);
var_opt.retrieveOption(argc,argv);
@@ -74,8 +83,10 @@ int main(int argc, char *argv[])
min_opt.retrieveOption(argc,argv);
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);
correlation_opt.retrieveOption(argc,argv);
rmse_opt.retrieveOption(argc,argv);
reg_opt.retrieveOption(argc,argv);
@@ -90,9 +101,20 @@ int main(int argc, char *argv[])
exit(0);//help was invoked, stop processing
}
+ statfactory::StatFactory stat;
+ if(rand_opt[0]>0){
+ gsl_rng* r=stat.getRandomGenerator(time(NULL));
+ //todo: init random number generator using time...
+ if(verbose_opt[0])
+ std::cout << "generating " << rand_opt[0] << " random numbers: " << std::endl;
+ for(unsigned int i=0;i<rand_opt[0];++i)
+ std::cout << i << " " << stat.getRandomValue(r,randdist_opt[0],randa_opt[0],randb_opt[0]) << std::endl;
+ }
vector< vector<double> > dataVector(col_opt.size());
- vector< vector<int> > statVector(col_opt.size());
+ vector< vector<double> > statVector(col_opt.size());
+ if(!input_opt.size())
+ exit(0);
FileReaderAscii asciiReader(input_opt[0]);
asciiReader.setFieldSeparator(fs_opt[0]);
asciiReader.setComment(comment_opt[0]);
@@ -101,7 +123,6 @@ int main(int argc, char *argv[])
asciiReader.setMaxRow(range_opt[1]);
asciiReader.readData(dataVector,col_opt);
assert(dataVector.size());
- statfactory::StatFactory stat;
double minValue=min_opt[0];
double maxValue=max_opt[0];
if(histogram_opt[0]){
@@ -112,6 +133,25 @@ int main(int argc, char *argv[])
nbin_opt[0]=maxValue-minValue+1;
}
}
+ double minX=min_opt[0];
+ double minY=(min_opt.size()==2)? min_opt[1] : min_opt[0];
+ double maxX=max_opt[0];
+ double maxY=(max_opt.size()==2)? max_opt[1] : max_opt[0];
+ if(histogram2d_opt[0]){
+ assert(col_opt.size()==2);
+ if(nbin_opt[0]<1){
+ std::cerr << "Warning: number of bins not defined, calculating bins from min and max value" << std::endl;
+ if(maxValue<=minValue){
+ stat.minmax(dataVector[0],dataVector[0].begin(),dataVector[0].end(),minX,maxX);
+ stat.minmax(dataVector[1],dataVector[1].begin(),dataVector[1].end(),minY,maxY);
+ }
+ minValue=(minX<minY)? minX:minY;
+ maxValue=(maxX>maxY)? maxX:maxY;
+ if(verbose_opt[0])
+ std::cout << "min and max values: " << minValue << ", " << maxValue << std::endl;
+ nbin_opt[0]=maxValue-minValue+1;
+ }
+ }
for(int icol=0;icol<col_opt.size();++icol){
if(!dataVector[icol].size()){
std::cerr << "Warning: dataVector[" << icol << "] is empty" << std::endl;
@@ -140,10 +180,22 @@ int main(int argc, char *argv[])
cout << "max value column " << col_opt[icol] << ": " << stat.max(dataVector[icol]) << endl;
}
if(histogram_opt[0]){
+ //todo: support kernel density function and estimate sigma as in practical estimate of the bandwith in http://en.wikipedia.org/wiki/Kernel_density_estimation
+ double sigma=0;
+ if(kde_opt.size()){
+ if(kde_opt[0]>0)
+ sigma=kde_opt[0];
+ else
+ sigma=1.06*sqrt(stat.var(dataVector[icol]))*pow(dataVector[icol].size(),-0.2);
+ }
assert(nbin_opt[0]);
- if(verbose_opt[0])
- std::cout << "calculating histogram for col " << icol << std::endl;
- stat.distribution(dataVector[icol],dataVector[icol].begin(),dataVector[icol].end(),statVector[icol],nbin_opt[0],minValue,maxValue);
+ if(verbose_opt[0]){
+ if(sigma>0)
+ std::cout << "calculating kernel density estimate with sigma " << sigma << " for col " << icol << std::endl;
+ else
+ std::cout << "calculating histogram for col " << icol << std::endl;
+ }
+ stat.distribution(dataVector[icol],dataVector[icol].begin(),dataVector[icol].end(),statVector[icol],nbin_opt[0],minValue,maxValue,sigma);
if(verbose_opt[0])
std::cout << "min and max values: " << minValue << ", " << maxValue << std::endl;
}
@@ -177,6 +229,38 @@ int main(int argc, char *argv[])
cout << endl;
}
}
+ if(histogram2d_opt[0]){
+ assert(nbin_opt[0]);
+ assert(dataVector.size()==2);
+ assert(dataVector[0].size()==dataVector[1].size());
+ double sigma=0;
+ //kernel density estimation as in http://en.wikipedia.org/wiki/Kernel_density_estimation
+ if(kde_opt.size()){
+ if(kde_opt[0]>0)
+ sigma=kde_opt[0];
+ else
+ sigma=1.06*sqrt(sqrt(stat.var(dataVector[0]))*sqrt(stat.var(dataVector[0])))*pow(dataVector[0].size(),-0.2);
+ }
+ assert(nbin_opt[0]);
+ if(verbose_opt[0]){
+ if(sigma>0)
+ std::cout << "calculating 2d kernel density estimate with sigma " << sigma << " for cols " << col_opt[0] << " and " << col_opt[1] << std::endl;
+ else
+ std::cout << "calculating 2d histogram for cols " << col_opt[0] << " and " << col_opt[1] << std::endl;
+ std::cout << "nbin: " << nbin_opt[0] << std::endl;
+ }
+ std::vector< std::vector<double> > histVector;
+ stat.distribution2d(dataVector[0],dataVector[1],histVector,nbin_opt[0],minX,maxX,minY,maxY,sigma);
+ for(int binX=0;binX<nbin_opt[0];++binX){
+ std::cout << std::endl;
+ for(int binY=0;binY<nbin_opt[0];++binY){
+ double value=0;
+ value=static_cast<double>(histVector[binX][binY])/dataVector[0].size();
+ std::cout << (maxX-minX)*binX/(nbin_opt[0]-1)+minX << " " << (maxY-minY)*binY/(nbin_opt[0]-1)+minY << " " << value << std::endl;
+ }
+ }
+ }
+
if(output_opt[0]){
for(int irow=0;irow<dataVector.begin()->size();++irow){
for(int icol=0;icol<col_opt.size();++icol){
@@ -187,4 +271,4 @@ int main(int argc, char *argv[])
cout << endl;
}
}
-}
+}
diff --git a/src/apps/pkstatogr.cc b/src/apps/pkstatogr.cc
index 871b2c5..81d8ec1 100644
--- a/src/apps/pkstatogr.cc
+++ b/src/apps/pkstatogr.cc
@@ -27,8 +27,8 @@ along with pktools. If not, see <http://www.gnu.org/licenses/>.
int main(int argc, char *argv[])
{
- Optionpk<string> input_opt("i", "input", "Input shape file", "");
- Optionpk<string> fieldname_opt("n", "fname", "fields on which to calculate statistics", "");
+ Optionpk<std::string> input_opt("i", "input", "Input shape file", "");
+ Optionpk<std::string> fieldname_opt("n", "fname", "fields on which to calculate statistics", "");
Optionpk<bool> minmax_opt("mm","minmax","calculate minimum and maximum value",false);
Optionpk<double> min_opt("min","min","set minimum value",0);
Optionpk<double> max_opt("max","max","set maximum value",0);
@@ -57,7 +57,7 @@ int main(int argc, char *argv[])
size_opt.retrieveOption(argc,argv);
verbose_opt.retrieveOption(argc,argv);
}
- catch(string predefinedString){
+ catch(std::string predefinedString){
std::cout << predefinedString << std::endl;
exit(0);
}
@@ -70,12 +70,12 @@ int main(int argc, char *argv[])
try{
imgReader.open(input_opt[0]);
}
- catch(string errorstring){
+ catch(std::string errorstring){
std::cerr << errorstring << std::endl;
}
ImgReaderOgr inputReader(input_opt[0]);
- vector<double> theData;
+ std::vector<double> theData;
statfactory::StatFactory stat;
//todo: implement ALL
@@ -84,7 +84,7 @@ int main(int argc, char *argv[])
std::cout << "field: " << ifield << std::endl;
theData.clear();
inputReader.readData(theData,OFTReal,fieldname_opt[ifield],0,verbose_opt[0]);
- vector<int> binData;
+ std::vector<double> binData;
double minValue=min_opt[0];
double maxValue=max_opt[0];
if(histogram_opt[0]){
@@ -97,7 +97,7 @@ int main(int argc, char *argv[])
try{
stat.distribution(theData,theData.begin(),theData.end(),binData,nbin_opt[0],minValue,maxValue);
}
- catch(string theError){
+ catch(std::string theError){
std::cerr << "Warning: all identical values in data" << std::endl;
exit(1);
}
@@ -113,8 +113,8 @@ int main(int argc, char *argv[])
if(stdev_opt[0])
std::cout << " --stdev " << sqrt(theVar);
if(minmax_opt[0]){
- cout << " -min " << stat.min(theData);
- cout << " -max " << stat.max(theData);
+ std::cout << " -min " << stat.min(theData);
+ std::cout << " -max " << stat.max(theData);
}
if(median_opt[0])
std::cout << " -median " << stat.median(theData);
@@ -130,15 +130,15 @@ int main(int argc, char *argv[])
}
}
}
- catch(string theError){
+ catch(std::string theError){
if(mean_opt[0])
std::cout << " --mean " << theData.back();
if(stdev_opt[0])
std::cout << " --stdev " << "0";
if(min_opt[0])
- cout << " -min " << theData.back();
+ std::cout << " -min " << theData.back();
if(max_opt[0])
- cout << " -max " << theData.back();
+ std::cout << " -max " << theData.back();
if(median_opt[0])
std::cout << " -median " << theData.back();
if(size_opt[0])
diff --git a/src/base/PosValue.h b/src/base/PosValue.h
index 3a9856f..86e08d8 100644
--- a/src/base/PosValue.h
+++ b/src/base/PosValue.h
@@ -20,8 +20,6 @@ along with pktools. If not, see <http://www.gnu.org/licenses/>.
#ifndef _POSVALUE_H_
#define _POSVALUE_H_
-using namespace std;
-
struct PosValue{
double posx;
double posy;
diff --git a/src/base/Vector2d.h b/src/base/Vector2d.h
index 3cb8703..a10c192 100644
--- a/src/base/Vector2d.h
+++ b/src/base/Vector2d.h
@@ -62,9 +62,9 @@ public:
void scale(const std::vector<double> &scaleVector, const std::vector<double> &offsetVector, Vector2d<T>& scaledOutput);
void scale(const T lbound, const T ubound, std::vector<double> &scaleVector, std::vector<double> &offsetVector, Vector2d<T>& scaledOutput);
Vector2d<T> operator=(const Vector2d<T>& v1);
-// ostream& operator<<(ostream& os, const Vector2d<T>& v);
-// template<class T> ostream& operator<<(ostream& os, const Vector2d<T>& v);
- template<class T1> friend ostream& operator<<(ostream & os, const Vector2d<T1>& v);
+// std::ostream& operator<<(std::ostream& os, const Vector2d<T>& v);
+// template<class T> std::ostream& operator<<(std::ostream& os, const Vector2d<T>& v);
+ template<class T1> friend std::ostream& operator<<(std::ostream & os, const Vector2d<T1>& v);
Vector2d<T> sum(const Vector2d<T>& v1, const Vector2d<T>& v2) const;
T max(int& x, int& y, double maxValue) const;
@@ -192,13 +192,13 @@ template<class T> void Vector2d<T>::selectCols(const std::list<int> &cols)
(*this)[irow].erase(((*this)[irow]).begin()+icol);
}
-template<class T1> ostream& operator<<(ostream& os, const Vector2d<T1>& v)
+template<class T1> std::ostream& operator<<(std::ostream& os, const Vector2d<T1>& v)
{
for(int irow=0;irow<v.size();++irow){
for(int icol=0;icol<v[irow].size();++icol){
os << v[irow][icol] << "\t";
}
- os << endl;
+ os << std::endl;
}
return os;
// os << theOption.getLongName() << ": ";
diff --git a/src/imageclasses/ImgReaderOgr.h b/src/imageclasses/ImgReaderOgr.h
index ecf06f3..27890fd 100644
--- a/src/imageclasses/ImgReaderOgr.h
+++ b/src/imageclasses/ImgReaderOgr.h
@@ -48,7 +48,7 @@ public:
template <typename T> int readData(Vector2d<T>& data, const OGRFieldType& fieldType, std::vector<std::string>& fields, int layer=0, bool pos=false, bool verbose=false);//default layer 0 and no pos information in data
template <typename T> int readData(std::map<int,Vector2d<T> >& data, const OGRFieldType& fieldType, std::vector<std::string>& fields, const std::string& label, int layer=0, bool pos=false, bool verbose=false);//default layer 0 and no pos information in data
template <typename T> int readData(std::map<std::string,Vector2d<T> >& data, const OGRFieldType& fieldType, std::vector<std::string>& fields, const std::string& label, int layer=0, bool pos=false, bool verbose=false);//default layer 0 and no pos information in data
- void shape2ascii(ostream& theOstream, const std::string& pointname, int layer=0, bool verbose=false);
+ void shape2ascii(std::ostream& theOstream, const std::string& pointname, int layer=0, bool verbose=false);
unsigned long int getFeatureCount(int layer=0) const;
int getFieldCount(int layer=0) const;
OGRLayer* getLayer(int layer=0){return m_datasource->GetLayer(layer);};
@@ -65,7 +65,7 @@ public:
template<typename T> int readSql(std::map<int,Vector2d<T> >& data, const OGRFieldType& fieldType, std::vector<std::string>& fields, const std::string& label, const std::string& sqlStatement, OGRGeometry* spatialFilter, int layer=0, bool pos=false, bool verbose=false);
bool getExtent(double& ulx, double& uly, double& lrx, double& lry, int layer=0);
- friend ostream& operator<<(ostream& theOstream, ImgReaderOgr& theImageReader);
+ friend std::ostream& operator<<(std::ostream& theOstream, ImgReaderOgr& theImageReader);
protected:
void setCodec(void);
@@ -82,36 +82,36 @@ template <typename T> int ImgReaderOgr::readData(std::map<int,Vector2d<T> >& dat
assert(m_datasource->GetLayerCount()>layer);
OGRLayer *poLayer;
if(verbose)
- cout << "number of layers: " << m_datasource->GetLayerCount() << endl;
+ std::cout << "number of layers: " << m_datasource->GetLayerCount() << std::endl;
poLayer = m_datasource->GetLayer(layer);
if(poLayer!=NULL){
OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
if(fields.empty()){
fields.resize(poFDefn->GetFieldCount());
if(verbose)
- cout << "resized fields to " << fields.size() << endl;
+ std::cout << "resized fields to " << fields.size() << std::endl;
}
//start reading features from the layer
OGRFeature *poFeature;
if(verbose)
- cout << "reset reading" << endl;
+ std::cout << "reset reading" << std::endl;
poLayer->ResetReading();
unsigned long int ifeature=0;
int posOffset=(pos)?2:0;
if(verbose)
- cout << "going through features" << endl << flush;
+ std::cout << "going through features" << std::endl << std::flush;
int theClass=0;
while( (poFeature = poLayer->GetNextFeature()) != NULL ){
std::vector<T> theFeature;//(fields.size()+posOffset);//x,y+selectedfields
if(verbose)
- cout << "reading feature " << ifeature << endl << flush;
+ std::cout << "reading feature " << ifeature << std::endl << std::flush;
OGRGeometry *poGeometry;
poGeometry = poFeature->GetGeometryRef();
if(verbose){
if(poGeometry == NULL)
- cerr << "no geometry defined" << endl << flush;
+ std::cerr << "no geometry defined" << std::endl << std::flush;
else if(wkbFlatten(poGeometry->getGeometryType()) != wkbPoint)
- cerr << "Warning: poGeometry type: " << wkbFlatten(poGeometry->getGeometryType()) << endl << flush;
+ std::cerr << "Warning: poGeometry type: " << wkbFlatten(poGeometry->getGeometryType()) << std::endl << std::flush;
}
assert(poGeometry != NULL );
// && wkbFlatten(poGeometry->getGeometryType()) == wkbPoint);
@@ -180,14 +180,14 @@ template <typename T> int ImgReaderOgr::readData(std::map<int,Vector2d<T> >& dat
++ifeature;
}
if(verbose)
- cout << "number of features read: " << ifeature << endl << flush;
+ std::cout << "number of features read: " << ifeature << std::endl << std::flush;
typename std::map<int,Vector2d<T> >::const_iterator mit=data.begin();
int nband=0;
if(verbose)
- cout << "read classes: " << flush;
+ std::cout << "read classes: " << std::flush;
while(mit!=data.end()){
if(verbose)
- cout << mit->first << " " << flush;
+ std::cout << mit->first << " " << std::flush;
if(!nband)
nband=fields.size();
if(pos)
@@ -197,7 +197,7 @@ template <typename T> int ImgReaderOgr::readData(std::map<int,Vector2d<T> >& dat
++mit;
}
if(verbose)
- cout << endl << flush;
+ std::cout << std::endl << std::flush;
return(nband);
}
else{
@@ -215,7 +215,7 @@ template <typename T> int ImgReaderOgr::readData(std::map<std::string,Vector2d<T
assert(m_datasource->GetLayerCount()>layer);
OGRLayer *poLayer;
if(verbose)
- cout << "number of layers: " << m_datasource->GetLayerCount() << endl;
+ std::cout << "number of layers: " << m_datasource->GetLayerCount() << std::endl;
poLayer = m_datasource->GetLayer(layer);
if(poLayer!=NULL){
OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
@@ -224,30 +224,30 @@ template <typename T> int ImgReaderOgr::readData(std::map<std::string,Vector2d<T
if(fields.empty()){
fields.resize(poFDefn->GetFieldCount());
if(verbose)
- cout << "resized fields to " << fields.size() << endl;
+ std::cout << "resized fields to " << fields.size() << std::endl;
}
//start reading features from the layer
OGRFeature *poFeature;
if(verbose)
- cout << "reset reading" << endl;
+ std::cout << "reset reading" << std::endl;
poLayer->ResetReading();
unsigned long int ifeature=0;
int posOffset=(pos)?2:0;
if(verbose)
- cout << "going through features to fill in string map" << endl << flush;
+ std::cout << "going through features to fill in string map" << std::endl << std::flush;
std::string theClass;
while( (poFeature = poLayer->GetNextFeature()) != NULL ){
std::vector<T> theFeature;//(fields.size()+posOffset);//x,y+selectedfields
if(verbose)
- cout << "reading feature " << ifeature << endl << flush;
+ std::cout << "reading feature " << ifeature << std::endl << std::flush;
OGRGeometry *poGeometry;
poGeometry = poFeature->GetGeometryRef();
if(verbose){
if(poGeometry == NULL)
- cerr << "no geometry defined" << endl << flush;
+ std::cerr << "no geometry defined" << std::endl << std::flush;
else if(wkbFlatten(poGeometry->getGeometryType()) != wkbPoint)
- cerr << "Warning: poGeometry type: " << wkbFlatten(poGeometry->getGeometryType()) << endl << flush;
+ std::cerr << "Warning: poGeometry type: " << wkbFlatten(poGeometry->getGeometryType()) << std::endl << std::flush;
}
assert(poGeometry != NULL );
// && wkbFlatten(poGeometry->getGeometryType()) == wkbPoint);
@@ -318,14 +318,14 @@ template <typename T> int ImgReaderOgr::readData(std::map<std::string,Vector2d<T
++ifeature;
}
if(verbose)
- cout << "number of features read: " << ifeature << endl << flush;
+ std::cout << "number of features read: " << ifeature << std::endl << std::flush;
typename std::map<std::string,Vector2d<T> >::const_iterator mit=data.begin();
int nband=0;
if(verbose)
- cout << "read classes: " << flush;
+ std::cout << "read classes: " << std::flush;
while(mit!=data.end()){
if(verbose)
- cout << mit->first << " " << flush;
+ std::cout << mit->first << " " << std::flush;
if(!nband)
nband=fields.size();
if(pos)
@@ -335,7 +335,7 @@ template <typename T> int ImgReaderOgr::readData(std::map<std::string,Vector2d<T
++mit;
}
if(verbose)
- cout << endl << flush;
+ std::cout << std::endl << std::flush;
return(nband);
}
else{
@@ -352,27 +352,27 @@ template <typename T> int ImgReaderOgr::readXY(std::vector<T>& xVector, std::vec
assert(m_datasource->GetLayerCount()>layer);
OGRLayer *poLayer;
if(verbose)
- cout << "number of layers: " << m_datasource->GetLayerCount() << endl;
+ std::cout << "number of layers: " << m_datasource->GetLayerCount() << std::endl;
poLayer = m_datasource->GetLayer(layer);
OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
//start reading features from the layer
OGRFeature *poFeature;
if(verbose)
- cout << "reset reading" << endl;
+ std::cout << "reset reading" << std::endl;
poLayer->ResetReading();
unsigned long int ifeature=0;
if(verbose)
- cout << "going through features" << endl << flush;
+ std::cout << "going through features" << std::endl << std::flush;
while( (poFeature = poLayer->GetNextFeature()) != NULL ){
if(verbose)
- cout << "reading feature " << ifeature << endl << flush;
+ std::cout << "reading feature " << ifeature << std::endl << std::flush;
OGRGeometry *poGeometry;
poGeometry = poFeature->GetGeometryRef();
if(verbose){
if(poGeometry == NULL)
- cerr << "no geometry defined" << endl << flush;
+ std::cerr << "no geometry defined" << std::endl << std::flush;
else// if(wkbFlatten(poGeometry->getGeometryType()) != wkbPoint)
- cout << "poGeometry type: " << wkbFlatten(poGeometry->getGeometryType()) << endl;
+ std::cout << "poGeometry type: " << wkbFlatten(poGeometry->getGeometryType()) << std::endl;
}
// assert(poGeometry != NULL
// && wkbFlatten(poGeometry->getGeometryType()) == wkbPoint);
@@ -400,21 +400,21 @@ template <typename T> int ImgReaderOgr::readData(std::vector<T>& data, const OGR
assert(m_datasource->GetLayerCount()>layer);
OGRLayer *poLayer;
if(verbose)
- cout << "number of layers: " << m_datasource->GetLayerCount() << endl;
+ std::cout << "number of layers: " << m_datasource->GetLayerCount() << std::endl;
poLayer = m_datasource->GetLayer(layer);
OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
if(fields.empty()){
fields.resize(poFDefn->GetFieldCount());
if(verbose)
- cout << "resized fields to " << fields.size() << endl;
+ std::cout << "resized fields to " << fields.size() << std::endl;
}
OGRGeometry *poGeometry;
poGeometry = poFeature->GetGeometryRef();
if(verbose){
if(poGeometry == NULL)
- cerr << "no geometry defined" << endl << flush;
+ std::cerr << "no geometry defined" << std::endl << std::flush;
else if(wkbFlatten(poGeometry->getGeometryType()) != wkbPoint)
- cerr << "poGeometry type: " << wkbFlatten(poGeometry->getGeometryType()) << endl << flush;
+ std::cerr << "poGeometry type: " << wkbFlatten(poGeometry->getGeometryType()) << std::endl << std::flush;
}
assert(poGeometry != NULL);
// && wkbFlatten(poGeometry->getGeometryType()) == wkbPoint);
@@ -497,36 +497,36 @@ template <typename T> int ImgReaderOgr::readData(std::vector<T>& data, const OGR
assert(m_datasource->GetLayerCount()>layer);
OGRLayer *poLayer;
if(verbose)
- cout << "number of layers: " << m_datasource->GetLayerCount() << endl;
+ std::cout << "number of layers: " << m_datasource->GetLayerCount() << std::endl;
poLayer = m_datasource->GetLayer(layer);
OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
int nfield=(theField!="")? poFDefn->GetFieldCount() : 1;
if(theField==""){
//read first field available
if(verbose)
- cout << "read first field from total of " << nfield << endl;
+ std::cout << "read first field from total of " << nfield << std::endl;
}
//start reading features from the layer
OGRFeature *poFeature;
if(verbose)
- cout << "reset reading" << endl;
+ std::cout << "reset reading" << std::endl;
poLayer->ResetReading();
unsigned long int ifeature=0;
if(verbose)
- cout << "going through features" << endl << flush;
+ std::cout << "going through features" << std::endl << std::flush;
while( (poFeature = poLayer->GetNextFeature()) != NULL ){
// std::vector<T> theFeature;//(fields.size()+posOffset);//x,y+selectedfields
T theFeature;
if(verbose)
- cout << "reading feature " << ifeature << endl << flush;
+ std::cout << "reading feature " << ifeature << std::endl << std::flush;
OGRGeometry *poGeometry;
poGeometry = poFeature->GetGeometryRef();
if(verbose){
if(poGeometry == NULL)
- cerr << "no geometry defined" << endl << flush;
+ std::cerr << "no geometry defined" << std::endl << std::flush;
else// if(wkbFlatten(poGeometry->getGeometryType()) != wkbPoint)
- cout << "poGeometry type: " << wkbFlatten(poGeometry->getGeometryType()) << endl;
+ std::cout << "poGeometry type: " << wkbFlatten(poGeometry->getGeometryType()) << std::endl;
}
// assert(poGeometry != NULL
// && wkbFlatten(poGeometry->getGeometryType()) == wkbPoint);
@@ -552,7 +552,7 @@ template <typename T> int ImgReaderOgr::readData(std::vector<T>& data, const OGR
}
data.push_back(theFeature);
if(verbose)
- cout << "feature is: " << theFeature << endl;
+ std::cout << "feature is: " << theFeature << std::endl;
++ifeature;
}
if(data.size()){
@@ -573,34 +573,34 @@ template <typename T> int ImgReaderOgr::readData(Vector2d<T>& data, const OGRFie
assert(m_datasource->GetLayerCount()>layer);
OGRLayer *poLayer;
if(verbose)
- cout << "number of layers: " << m_datasource->GetLayerCount() << endl;
+ std::cout << "number of layers: " << m_datasource->GetLayerCount() << std::endl;
poLayer = m_datasource->GetLayer(layer);
OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
if(fields.empty()){
fields.resize(poFDefn->GetFieldCount());
if(verbose)
- cout << "resized fields to " << fields.size() << endl;
+ std::cout << "resized fields to " << fields.size() << std::endl;
}
//start reading features from the layer
OGRFeature *poFeature;
if(verbose)
- cout << "reset reading" << endl;
+ std::cout << "reset reading" << std::endl;
poLayer->ResetReading();
unsigned long int ifeature=0;
int posOffset=(pos)?2:0;
if(verbose)
- cout << "going through features" << endl << flush;
+ std::cout << "going through features" << std::endl << std::flush;
while( (poFeature = poLayer->GetNextFeature()) != NULL ){
std::vector<T> theFeature;//(fields.size()+posOffset);//x,y+selectedfields
if(verbose)
- cout << "reading feature " << ifeature << endl << flush;
+ std::cout << "reading feature " << ifeature << std::endl << std::flush;
OGRGeometry *poGeometry;
poGeometry = poFeature->GetGeometryRef();
if(verbose){
if(poGeometry == NULL)
- cerr << "no geometry defined" << endl << flush;
+ std::cerr << "no geometry defined" << std::endl << std::flush;
else if(wkbFlatten(poGeometry->getGeometryType()) != wkbPoint)
- cerr << "poGeometry type: " << wkbFlatten(poGeometry->getGeometryType()) << endl << flush;
+ std::cerr << "poGeometry type: " << wkbFlatten(poGeometry->getGeometryType()) << std::endl << std::flush;
}
assert(poGeometry != NULL
&& wkbFlatten(poGeometry->getGeometryType()) == wkbPoint);
@@ -673,29 +673,29 @@ template<typename T> int ImgReaderOgr::readSql(std::map<int, Vector2d<T> >& data
if(fields.empty()){
fields.resize(poFDefn->GetFieldCount());
if(verbose)
- cout << "resized fields to " << fields.size() << endl;
+ std::cout << "resized fields to " << fields.size() << std::endl;
}
//start reading features from the layer
OGRFeature *poFeature;
if(verbose)
- cout << "reset reading" << endl;
+ std::cout << "reset reading" << std::endl;
poLayer->ResetReading();
unsigned long int ifeature=0;
int posOffset=(pos)?2:0;
if(verbose)
- cout << "going through features" << endl << flush;
+ std::cout << "going through features" << std::endl << std::flush;
int theClass=0;
while( (poFeature = poLayer->GetNextFeature()) != NULL ){
std::vector<T> theFeature;//(fields.size()+posOffset);//x,y+selectedfields
if(verbose)
- cout << "reading feature " << ifeature << endl << flush;
+ std::cout << "reading feature " << ifeature << std::endl << std::flush;
OGRGeometry *poGeometry;
poGeometry = poFeature->GetGeometryRef();
if(verbose){
if(poGeometry == NULL)
- cerr << "no geometry defined" << endl << flush;
+ std::cerr << "no geometry defined" << std::endl << std::flush;
else if(wkbFlatten(poGeometry->getGeometryType()) != wkbPoint)
- cerr << "poGeometry type: " << wkbFlatten(poGeometry->getGeometryType()) << endl << flush;
+ std::cerr << "poGeometry type: " << wkbFlatten(poGeometry->getGeometryType()) << std::endl << std::flush;
}
assert(poGeometry != NULL
&& wkbFlatten(poGeometry->getGeometryType()) == wkbPoint);
@@ -746,14 +746,14 @@ template<typename T> int ImgReaderOgr::readSql(std::map<int, Vector2d<T> >& data
++ifeature;
}
if(verbose)
- cout << "number of features read: " << ifeature << endl << flush;
+ std::cout << "number of features read: " << ifeature << std::endl << std::flush;
typename std::map<int,Vector2d<T> >::const_iterator mit=data.begin();
int nband=0;
if(verbose)
- cout << "read classes: " << flush;
+ std::cout << "read classes: " << std::flush;
while(mit!=data.end()){
if(verbose)
- cout << mit->first << " " << flush;
+ std::cout << mit->first << " " << std::flush;
if(!nband)
nband=fields.size();
if(pos)
@@ -763,7 +763,7 @@ template<typename T> int ImgReaderOgr::readSql(std::map<int, Vector2d<T> >& data
++mit;
}
if(verbose)
- cout << endl << flush;
+ std::cout << std::endl << std::flush;
return(nband);
}
else{
@@ -785,28 +785,28 @@ template<typename T> int ImgReaderOgr::readSql(Vector2d<T>& data, const OGRField
if(fields.empty()){
fields.resize(poFDefn->GetFieldCount());
if(verbose)
- cout << "resized fields to " << fields.size() << endl;
+ std::cout << "resized fields to " << fields.size() << std::endl;
}
//start reading features from the layer
OGRFeature *poFeature;
if(verbose)
- cout << "reset reading" << endl;
+ std::cout << "reset reading" << std::endl;
poLayer->ResetReading();
unsigned long int ifeature=0;
int posOffset=(pos)?2:0;
if(verbose)
- cout << "going through features" << endl << flush;
+ std::cout << "going through features" << std::endl << std::flush;
while( (poFeature = poLayer->GetNextFeature()) != NULL ){
std::vector<T> theFeature;//(fields.size()+posOffset);//x,y+selectedfields
if(verbose)
- cout << "reading feature " << ifeature << endl << flush;
+ std::cout << "reading feature " << ifeature << std::endl << std::flush;
OGRGeometry *poGeometry;
poGeometry = poFeature->GetGeometryRef();
if(verbose){
if(poGeometry == NULL)
- cerr << "no geometry defined" << endl << flush;
+ std::cerr << "no geometry defined" << std::endl << std::flush;
else if(wkbFlatten(poGeometry->getGeometryType()) != wkbPoint)
- cerr << "poGeometry type: " << wkbFlatten(poGeometry->getGeometryType()) << endl << flush;
+ std::cerr << "poGeometry type: " << wkbFlatten(poGeometry->getGeometryType()) << std::endl << std::flush;
}
assert(poGeometry != NULL
&& wkbFlatten(poGeometry->getGeometryType()) == wkbPoint);
--
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