[pktools] 366/375: nodata values for spectral/temporal filtering
Bas Couwenberg
sebastic at xs4all.nl
Wed Dec 3 21:54:30 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 954843e7916289eb57fc107234a1a3e529446b01
Author: Pieter Kempeneers <kempenep at gmail.com>
Date: Thu Nov 27 18:03:23 2014 +0100
nodata values for spectral/temporal filtering
---
ChangeLog | 7 +
configure.ac | 4 +-
pktools.pc | 2 +-
src/algorithms/Filter.cc | 8 ++
src/algorithms/Filter.h | 11 +-
src/algorithms/StatFactory.h | 336 +++++++++++++++++++++++++++++++------------
src/apps/pkfilter.cc | 1 +
7 files changed, 273 insertions(+), 96 deletions(-)
diff --git a/ChangeLog b/ChangeLog
index 5b9e617..d73fa64 100755
--- a/ChangeLog
+++ b/ChangeLog
@@ -322,6 +322,13 @@ version 2.5.4
- ImgWriteOgr
overwrite existing ogr datasets per default
+version 2.6.1
+ - API
+ Filter.h support nodata values
+ - pkfilter
+ Declare nodata option as double (see ticket #43500)
+ Support nodata values for filtering in spectral/temporal domain (see ticket #43713)
+
Todo:
- todo for API
ImgReaderGdal (ImgWriterGdal) open in update mode (check gdal_edit.py: http://searchcode.com/codesearch/view/18938404)
diff --git a/configure.ac b/configure.ac
index 8ef471e..469ba42 100644
--- a/configure.ac
+++ b/configure.ac
@@ -1,4 +1,4 @@
-AC_INIT([pktools], [2.5.4], [kempenep at gmail.com])
+AC_INIT([pktools], [2.6.1], [kempenep at gmail.com])
#AM_INIT_AUTOMAKE([-Wall -Werror foreign])
AM_INIT_AUTOMAKE([-Wall -Wno-extra-portability foreign])
#AM_INIT_AUTOMAKE([subdir-objects]) #not working due to bug in autoconf, see Debian list: Bug #752993)
@@ -97,7 +97,7 @@ AC_SUBST([LIBS])
# For information on how to properly maintain the library version information,
# refer to the libtool manual, section "Updating library version information":
# http://www.gnu.org/software/libtool/manual/html_node/Updating-version-info.html
-AC_SUBST([PKTOOLS_SO_VERSION], [1:2:0])
+AC_SUBST([PKTOOLS_SO_VERSION], [1:3:0])
# files to generate via autotools (.am or .in source files)
AC_CONFIG_HEADERS([config.h])
diff --git a/pktools.pc b/pktools.pc
index 996f47d..dc94e8c 100644
--- a/pktools.pc
+++ b/pktools.pc
@@ -6,6 +6,6 @@ includedir=${prefix}/include
Name: pktools
Description: API library for pktools
Requires: gdal gsl
-Version: 2.5.4
+Version: 2.6.1
Libs: -L${libdir} -lbase -lalgorithms -limageClasses -lfileClasses -llasClasses
Cflags: -I${includedir}/pktools
diff --git a/src/algorithms/Filter.cc b/src/algorithms/Filter.cc
index 4241d41..f2b6aa0 100644
--- a/src/algorithms/Filter.cc
+++ b/src/algorithms/Filter.cc
@@ -51,6 +51,13 @@ void filter::Filter::setTaps(const vector<double> &taps, bool normalize)
assert(m_taps.size()%2);
}
+int filter::Filter::pushNoDataValue(double noDataValue)
+{
+ if(find(m_noDataValues.begin(),m_noDataValues.end(),noDataValue)==m_noDataValues.end())
+ m_noDataValues.push_back(noDataValue);
+ return(m_noDataValues.size());
+}
+
void filter::Filter::dwtForward(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& wavelet_type, int family){
const char* pszMessage;
void* pProgressArg=NULL;
@@ -338,6 +345,7 @@ void filter::Filter::stat(const ImgReaderGdal& input, ImgWriterGdal& output, con
assert(output.nrOfCol()==input.nrOfCol());
vector<double> lineOutput(output.nrOfCol());
statfactory::StatFactory stat;
+ stat.setNoDataValues(m_noDataValues);
const char* pszMessage;
void* pProgressArg=NULL;
GDALProgressFunc pfnProgress=GDALTermProgress;
diff --git a/src/algorithms/Filter.h b/src/algorithms/Filter.h
index 705bfc6..a75626c 100644
--- a/src/algorithms/Filter.h
+++ b/src/algorithms/Filter.h
@@ -65,6 +65,7 @@ public:
void setTaps(const std::vector<double> &taps, bool normalize=true);
void pushClass(short theClass=1){m_class.push_back(theClass);};
void pushMask(short theMask=0){m_mask.push_back(theMask);};
+ int pushNoDataValue(double noDataValue=0);//{m_mask.push_back(theMask);};
template<class T> void filter(const std::vector<T>& input, std::vector<T>& output);
template<class T> void filter(const std::vector<T>& input, std::vector<T>& output, const std::string& method, int dim);
template<class T> void smooth(const std::vector<T>& input, std::vector<T>& output, short dim);
@@ -140,7 +141,8 @@ private:
std::vector<double> m_taps;
std::vector<short> m_class;
std::vector<short> m_mask;
- std::string m_padding;
+ std::string m_padding;
+ std::vector<double> m_noDataValues;
};
@@ -424,7 +426,7 @@ template<class T> void Filter::filter(const std::vector<T>& input, std::vector<T
int i=0;
//start: extend input by padding
for(i=0;i<m_taps.size()/2;++i){
- //todo:introduce nodata
+ //todo:introduce nodata?
output[i]=m_taps[m_taps.size()/2]*input[i];
for(int t=1;t<=m_taps.size()/2;++t){
output[i]+=m_taps[m_taps.size()/2+t]*input[i+t];
@@ -460,9 +462,9 @@ template<class T> void Filter::filter(const std::vector<T>& input, std::vector<T
}
//end: extend input by padding
for(i=input.size()-m_taps.size()/2;i<input.size();++i){
- //todo:introduce nodata
+ //todo:introduce nodata?
output[i]=m_taps[m_taps.size()/2]*input[i];
- //todo:introduce nodata
+ //todo:introduce nodata?
for(int t=1;t<=m_taps.size()/2;++t){
output[i]+=m_taps[m_taps.size()/2-t]*input[i-t];
if(i+t<input.size())
@@ -497,6 +499,7 @@ template<class T> void Filter::filter(const std::vector<T>& input, std::vector<T
output.resize(input.size());
int i=0;
statfactory::StatFactory stat;
+ stat.setNoDataValues(m_noDataValues);
std::vector<T> statBuffer;
short binValue=0;
//start: extend input by padding
diff --git a/src/algorithms/StatFactory.h b/src/algorithms/StatFactory.h
index abf7640..4e575e3 100644
--- a/src/algorithms/StatFactory.h
+++ b/src/algorithms/StatFactory.h
@@ -210,149 +210,238 @@ private:
template<class T> inline typename std::vector<T>::const_iterator StatFactory::mymin(const std::vector<T>& v, typename std::vector<T>::const_iterator begin, typename std::vector<T>::const_iterator end) const
{
+ bool isValid=false;
typename std::vector<T>::const_iterator tmpIt=begin;
- for (typename std::vector<T>::const_iterator it = begin; it!=end; ++it){
- if(!isNoData(*it))
+ for(typename std::vector<T>::const_iterator it = begin; it!=end; ++it){
+ if(!isNoData(*it)){
+ isValid=true;
if(*tmpIt>*it)
tmpIt=it;
+ }
+ }
+ if(isValid)
+ return tmpIt;
+ else if(m_noDataValues.size())
+ return m_noDataValues[0];
+ else{
+ std::string errorString="Error: no valid data found";
+ throw(errorString);
}
- return tmpIt;
}
template<class T> inline typename std::vector<T>::iterator StatFactory::mymin(const std::vector<T>& v, typename std::vector<T>::iterator begin, typename std::vector<T>::iterator end) const
{
+ bool isValid=false;
typename std::vector<T>::iterator tmpIt=begin;
for (typename std::vector<T>::const_iterator it = begin; it!=end; ++it){
- if(!isNoData(*it))
+ if(!isNoData(*it)){
+ isValid=true;
if(*tmpIt>*it)
tmpIt=it;
+ }
+ }
+ if(isValid)
+ return tmpIt;
+ else if(m_noDataValues.size())
+ return m_noDataValues[0];
+ else{
+ std::string errorString="Error: no valid data found";
+ throw(errorString);
}
- return tmpIt;
}
template<class T> inline typename std::vector<T>::const_iterator StatFactory::mymin(const std::vector<T>& v, typename std::vector<T>::const_iterator begin, typename std::vector<T>::const_iterator end, T minConstraint) const
{
+ bool isValid=false;
typename std::vector<T>::const_iterator tmpIt=v.end();
T minValue=minConstraint;
for (typename std::vector<T>::const_iterator it = begin; it!=end; ++it){
if(isNoData(*it))
continue;
+ isValid=true;
if((minConstraint<=*it)&&(*it<=minValue)){
tmpIt=it;
minValue=*it;
}
}
- return tmpIt;
+ if(isValid)
+ return tmpIt;
+ else if(m_noDataValues.size())
+ return m_noDataValues[0];
+ else{
+ std::string errorString="Error: no valid data found";
+ throw(errorString);
+ }
}
template<class T> inline typename std::vector<T>::iterator StatFactory::mymin(const std::vector<T>& v, typename std::vector<T>::iterator begin, typename std::vector<T>::iterator end, T minConstraint) const
{
+ bool isValid=false;
typename std::vector<T>::iterator tmpIt=v.end();
T minValue=minConstraint;
for (typename std::vector<T>::const_iterator it = begin; it!=end; ++it){
if(isNoData(*it))
continue;
+ isValid=true;
if((minConstraint<=*it)&&(*it<=minValue)){
tmpIt=it;
minValue=*it;
}
}
- return tmpIt;
+ if(isValid)
+ return tmpIt;
+ else if(m_noDataValues.size())
+ return m_noDataValues[0];
+ else{
+ std::string errorString="Error: no valid data found";
+ throw(errorString);
+ }
}
template<class T> inline typename std::vector<T>::const_iterator StatFactory::mymax(const std::vector<T>& v, typename std::vector<T>::const_iterator begin, typename std::vector<T>::const_iterator end) const
{
+ bool isValid=false;
typename std::vector<T>::const_iterator tmpIt=begin;
for (typename std::vector<T>::iterator it = begin; it!=end; ++it){
if(isNoData(*it))
continue;
+ isValid=true;
if(*tmpIt<*it)
tmpIt=it;
}
- return tmpIt;
+ if(isValid)
+ return tmpIt;
+ else if(m_noDataValues.size())
+ return m_noDataValues[0];
+ else{
+ std::string errorString="Error: no valid data found";
+ throw(errorString);
+ }
}
template<class T> inline typename std::vector<T>::iterator StatFactory::mymax(const std::vector<T>& v, typename std::vector<T>::iterator begin, typename std::vector<T>::iterator end) const
{
+ bool isValid=false;
typename std::vector<T>::iterator tmpIt=begin;
for (typename std::vector<T>::iterator it = begin; it!=end; ++it){
if(isNoData(*it))
continue;
+ isValid=true;
if(*tmpIt<*it)
tmpIt=it;
}
- return tmpIt;
+ if(isValid)
+ return tmpIt;
+ else
+ return end;
}
template<class T> inline typename std::vector<T>::const_iterator StatFactory::mymax(const std::vector<T>& v, typename std::vector<T>::const_iterator begin, typename std::vector<T>::const_iterator end, T maxConstraint) const
{
+ bool isValid=false;
typename std::vector<T>::const_iterator tmpIt=v.end();
T maxValue=maxConstraint;
for (typename std::vector<T>::const_iterator it = begin; it!=end; ++it){
if(isNoData(*it))
continue;
+ isValid=true;
if((maxValue<=*it)&&(*it<=maxConstraint)){
tmpIt=it;
maxValue=*it;
}
}
- return tmpIt;
+ if(isValid)
+ return tmpIt;
+ else
+ return end;
}
template<class T> inline typename std::vector<T>::iterator StatFactory::mymax(const std::vector<T>& v, typename std::vector<T>::iterator begin, typename std::vector<T>::iterator end, T maxConstraint) const
{
+ bool isValid=false;
typename std::vector<T>::iterator tmpIt=v.end();
T maxValue=maxConstraint;
for (typename std::vector<T>::const_iterator it = begin; it!=end; ++it){
if(isNoData(*it))
continue;
+ isValid=true;
if((maxValue<=*it)&&(*it<=maxConstraint)){
tmpIt=it;
maxValue=*it;
}
}
- return tmpIt;
+ if(isValid)
+ return tmpIt;
+ else
+ return end;
}
-
-
-
template<class T> inline T StatFactory::mymin(const std::vector<T>& v) const
{
+ bool isValid=false;
T minValue=*(v.begin());
for (typename std::vector<T>::const_iterator it = v.begin(); it!=v.end(); ++it){
if(isNoData(*it))
continue;
+ isValid=true;
if(minValue>*it)
minValue=*it;
}
- return minValue;
+ if(isValid)
+ return minValue;
+ else if(m_noDataValues.size())
+ return m_noDataValues[0];
+ else{
+ std::string errorString="Error: no valid data found";
+ throw(errorString);
+ }
}
template<class T> inline T StatFactory::mymin(const std::vector<T>& v, T minConstraint) const
{
+ bool isValid=false;
T minValue=minConstraint;
for (typename std::vector<T>::const_iterator it = v.begin(); it!=v.end(); ++it){
+ if(isNoData(*it))
+ continue;
+ isValid=true;
if((minConstraint<=*it)&&(*it<=minValue))
minValue=*it;
}
- return minValue;
+ if(isValid)
+ return minValue;
+ else if(m_noDataValues.size())
+ return m_noDataValues[0];
+ else{
+ std::string errorString="Error: no valid data found";
+ throw(errorString);
+ }
}
template<class T> inline T StatFactory::mymax(const std::vector<T>& v) const
{
+ bool isValid=false;
T maxValue=*(v.begin());
for (typename std::vector<T>::const_iterator it = v.begin(); it!=v.end(); ++it){
if(isNoData(*it))
continue;
+ isValid=true;
if(maxValue<*it)
maxValue=*it;
}
- return maxValue;
+ if(isValid)
+ return maxValue;
+ else if(m_noDataValues.size())
+ return m_noDataValues[0];
+ else{
+ std::string errorString="Error: no valid data found";
+ throw(errorString);
+ }
}
template<class T> inline T StatFactory::mymax(const std::vector<T>& v, T maxConstraint) const
{
+ bool isValid=false;
T maxValue=maxConstraint;
for (typename std::vector<T>::const_iterator it = v.begin(); it!=v.end(); ++it){
if(isNoData(*it))
@@ -360,56 +449,95 @@ template<class T> inline T StatFactory::mymax(const std::vector<T>& v, T maxCons
if((maxValue<=*it)&&(*it<=maxConstraint))
maxValue=*it;
}
- return maxValue;
+ if(isValid)
+ return maxValue;
+ else if(m_noDataValues.size())
+ return m_noDataValues[0];
+ else{
+ std::string errorString="Error: no valid data found";
+ throw(errorString);
+ }
}
template<class T> inline 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
{
+ bool isValid=false;
typename std::vector<T>::const_iterator tmpIt=begin;
for (typename std::vector<T>::const_iterator it = begin; it!=end; ++it){
if(isNoData(*it))
continue;
+ isValid=true;
if(abs(*tmpIt)<abs(*it))
tmpIt=it;
}
- return tmpIt;
+ if(isValid)
+ return tmpIt;
+ else
+ return end;
}
template<class T> inline 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
{
+ bool isValid=false;
typename std::vector<T>::const_iterator tmpIt=begin;
for (typename std::vector<T>::const_iterator it = begin; it!=end; ++it){
if(isNoData(*it))
continue;
+ isValid=true;
if(abs(*tmpIt)>abs(*it))
tmpIt=it;
}
+ if(isValid)
+ return tmpIt;
+ else
+ return end;
}
template<class T> inline 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
{
+ bool isValid=false;
theMin=*begin;
theMax=*begin;
for (typename std::vector<T>::const_iterator it = begin; it!=end; ++it){
if(isNoData(*it))
continue;
+ isValid=true;
if(theMin>*it)
theMin=*it;
if(theMax<*it)
theMax=*it;
}
+ if(!isValid){
+ if(m_noDataValues.size()){
+ theMin=m_noDataValues[0];
+ theMax=m_noDataValues[0];
+ }
+ else{
+ std::string errorString="Error: no valid data found";
+ throw(errorString);
+ }
+ }
}
template<class T> inline T StatFactory::sum(const std::vector<T>& v) const
{
+ bool isValid=false;
typename std::vector<T>::const_iterator it;
T tmpSum=0;
for (it = v.begin(); it!= v.end(); ++it){
if(isNoData(*it))
continue;
+ isValid=true;
tmpSum+=*it;
}
- return tmpSum;
+ if(isValid)
+ return tmpSum;
+ else if(m_noDataValues.size())
+ return m_noDataValues[0];
+ else{
+ std::string errorString="Error: no valid data found";
+ throw(errorString);
+ }
}
template<class T> inline double StatFactory::mean(const std::vector<T>& v) const
@@ -426,18 +554,24 @@ template<class T> inline double StatFactory::mean(const std::vector<T>& v) const
}
if(validSize)
return static_cast<double>(tmpSum)/validSize;
- else
- return 0;
+ else if(m_noDataValues.size())
+ return m_noDataValues[0];
+ else{
+ std::string errorString="Error: no valid data found";
+ throw(errorString);
+ }
}
template<class T> inline void StatFactory::eraseNoData(std::vector<T>& v) const
{
- typename std::vector<T>::iterator it=v.begin();
- while(it!=v.end()){
- if(isNoData(*it))
- v.erase(it);
- else
- ++it;
+ if(m_noDataValues.size()){
+ typename std::vector<T>::iterator it=v.begin();
+ while(it!=v.end()){
+ if(isNoData(*it))
+ v.erase(it);
+ else
+ ++it;
+ }
}
}
@@ -445,11 +579,19 @@ template<class T> T StatFactory::median(const std::vector<T>& v) const
{
std::vector<T> tmpV=v;
eraseNoData(tmpV);
- sort(tmpV.begin(),tmpV.end());
- if(tmpV.size()%2)
- return tmpV[tmpV.size()/2];
- else
- return 0.5*(tmpV[tmpV.size()/2-1]+tmpV[tmpV.size()/2]);
+ if(tmpV.size()){
+ sort(tmpV.begin(),tmpV.end());
+ if(tmpV.size()%2)
+ return tmpV[tmpV.size()/2];
+ else
+ return 0.5*(tmpV[tmpV.size()/2-1]+tmpV[tmpV.size()/2]);
+ }
+ else if(m_noDataValues.size())
+ return m_noDataValues[0];
+ else{
+ std::string errorString="Error: no valid data found";
+ throw(errorString);
+ }
}
template<class T> double StatFactory::var(const std::vector<T>& v) const
@@ -470,17 +612,12 @@ template<class T> double StatFactory::var(const std::vector<T>& v) const
m1/=validSize;
return m2-m1*m1;
}
- else
- return 0;
- /* double v1=0; */
- /* double m1=mean(v); */
- /* double n=v.size(); */
- /* assert(n>1); */
- /* for (it = v.begin(); it!= v.end(); ++it) */
- /* v1+=(*it-m1)*(*it-m1); */
- /* v1/=(n-1); */
- /* assert(v1>=0); */
- /* return v1; */
+ else if(m_noDataValues.size())
+ return m_noDataValues[0];
+ else{
+ std::string errorString="Error: no valid data found";
+ throw(errorString);
+ }
}
template<class T> double StatFactory::moment(const std::vector<T>& v, int n) const
@@ -498,9 +635,12 @@ template<class T> double StatFactory::moment(const std::vector<T>& v, int n) con
}
if(validSize)
return m/validSize;
- else
- return 0;
- /* return m/v.size(); */
+ else if(m_noDataValues.size())
+ return m_noDataValues[0];
+ else{
+ std::string errorString="Error: no valid data found";
+ throw(errorString);
+ }
}
//central moment
@@ -517,17 +657,25 @@ template<class T> double StatFactory::cmoment(const std::vector<T>& v, int n) co
m+=pow((*it-m1),n);
++validSize;
}
- return m/validSize;
- return m/v.size();
+ if(validSize)
+ return m/validSize;
+ else if(m_noDataValues.size())
+ return m_noDataValues[0];
+ else{
+ std::string errorString="Error: no valid data found";
+ throw(errorString);
+ }
}
template<class T> double StatFactory::skewness(const std::vector<T>& v) const
{
+ //todo: what if nodata value?
return cmoment(v,3)/pow(var(v),1.5);
}
template<class T> double StatFactory::kurtosis(const std::vector<T>& v) const
{
+ //todo: what if nodata value?
double m2=cmoment(v,2);
double m4=cmoment(v,4);
return m4/m2/m2-3.0;
@@ -552,15 +700,14 @@ template<class T> void StatFactory::meanVar(const std::vector<T>& v, double& m1,
m1/=validSize;
v1=m2-m1*m1;
}
- /* typename std::vector<T>::const_iterator it; */
- /* v1=0; */
- /* m1=mean(v); */
- /* double n=v.size(); */
- /* assert(n>1); */
- /* for (it = v.begin(); it!= v.end(); ++it) */
- /* v1+=(*(it)-m1)*(*(it)-m1); */
- /* v1/=(n-1); */
- /* assert(v1>=0); */
+ else if(m_noDataValues.size()){
+ m1=m_noDataValues[0];
+ v1=m_noDataValues[0];
+ }
+ else{
+ std::string errorString="Error: no valid data found";
+ throw(errorString);
+ }
}
template<class T1, class T2> void StatFactory::scale2byte(const std::vector<T1>& input, std::vector<T2>& output, unsigned char lbound, unsigned char ubound) const
@@ -570,8 +717,10 @@ template<class T1, class T2> void StatFactory::scale2byte(const std::vector<T1>&
T1 maximum=mymax(input);
assert(maximum>minimum);
double scale=(ubound-lbound)/(maximum-minimum);
- for (int i=0;i<input.size();++i)
+ //todo: what if nodata value?
+ for (int i=0;i<input.size();++i){
output[i]=scale*(input[i]-(minimum))+lbound;
+ }
}
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) const
@@ -594,11 +743,15 @@ template<class T> void StatFactory::distribution(const std::vector<T>& input, t
throw(s.str());
}
assert(nbin);
- assert(input.size());
+ if(!input.size()){
+ std::string errorString="Error: no valid data found";
+ throw(errorString);
+ }
if(output.size()!=nbin){
output.resize(nbin);
for(int i=0;i<nbin;output[i++]=0);
}
+ bool isValid=false;
typename std::vector<T>::const_iterator it;
for(it=begin;it!=end;++it){
if(*it<minimum)
@@ -607,6 +760,7 @@ template<class T> void StatFactory::distribution(const std::vector<T>& input, t
continue;
if(isNoData(*it))
continue;
+ isValid=true;
if(sigma>0){
// minimum-=2*sigma;
// maximum+=2*sigma;
@@ -632,7 +786,11 @@ template<class T> void StatFactory::distribution(const std::vector<T>& input, t
// ++output[static_cast<int>(static_cast<double>((*it)-minimum)/(maximum-minimum)*nbin)];
}
}
- if(!filename.empty()){
+ if(!isValid){
+ std::string errorString="Error: no valid data found";
+ throw(errorString);
+ }
+ else if(!filename.empty()){
std::ofstream outputfile;
outputfile.open(filename.c_str());
if(!outputfile){
@@ -726,14 +884,11 @@ template<class T> void StatFactory::distribution2d(const std::vector<T>& inputX
}
}
+//todo: what with nodata values?
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) const
{
if(maximum<=minimum)
minmax(input,begin,end,minimum,maximum);
- // if(!minimum)
- // minimum=*(min(input,begin,end));
- // if(!maximum)
- // maximum=*(max(input,begin,end));
assert(maximum>minimum);
assert(nbin>1);
assert(input.size());
@@ -757,7 +912,6 @@ template<class T> void StatFactory::percentiles (const std::vector<T>& input, t
++vit;
}
if(inputBin.size()){
- /* output[ibin]=median(inputBin); */
output[ibin]=mymax(inputBin);
}
}
@@ -775,27 +929,6 @@ template<class T> void StatFactory::percentiles (const std::vector<T>& input, t
}
}
-// 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(std::vector<int>::iterator it=begin+1;it!=end;++it)
-// *it+=*(it-1);
-// if(!filename.empty()){
-// std::ofstream outputfile;
-// outputfile.open(filename.c_str());
-// if(!outputfile){
-// 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() << std::endl;
-// outputfile.close();
-// }
-// }
-
template<class T> void StatFactory::signature(const std::vector<T>& input, double&k, double& alpha, double& beta, double e) const
{
double m1=moment(input,1);
@@ -803,6 +936,7 @@ template<class T> void StatFactory::signature(const std::vector<T>& input, doubl
signature(m1,m2,k,alpha,beta,e);
}
+//todo: what with nodata values?
template<class T> void StatFactory::normalize(const std::vector<T>& input, std::vector<double>& output) const{
double total=sum(input);
if(total){
@@ -814,6 +948,7 @@ template<class T> void StatFactory::normalize(const std::vector<T>& input, std::
output=input;
}
+//todo: what with nodata values?
template<class T> void StatFactory::normalize_pct(std::vector<T>& input) const{
double total=sum(input);
if(total){
@@ -828,6 +963,8 @@ template<class T> double StatFactory::rmse(const std::vector<T>& x, const std::v
assert(x.size());
double mse=0;
for(int isample=0;isample<x.size();++isample){
+ if(isNoData(x[isample])||isNoData(y[isample]))
+ continue;
double e=x[isample]-y[isample];
mse+=e*e/x.size();
}
@@ -843,6 +980,7 @@ template<class T> double StatFactory::correlation(const std::vector<T>& x, const
meanVar(x,meanX,varX);
meanVar(y,meanY,varY);
double denom = sqrt(varX*varY);
+ bool isValid=false;
if(denom){
//Calculate the correlation series
sXY = 0;
@@ -850,19 +988,31 @@ template<class T> double StatFactory::correlation(const std::vector<T>& x, const
int j = i + delay;
if (j < 0 || j >= y.size())
continue;
+ else if(isNoData(x[i])||isNoData(y[j]))
+ continue;
else{
+ isValid=true;
assert(i>=0&&i<x.size());
assert(j>=0&&j<y.size());
sXY += (x[i] - meanX) * (y[j] - meanY);
}
}
- double minSize=(x.size()<y.size())?x.size():y.size();
- return(sXY / denom / (minSize-1));
+ if(isValid){
+ double minSize=(x.size()<y.size())?x.size():y.size();
+ return(sXY / denom / (minSize-1));
+ }
+ else if(m_noDataValues.size())
+ return m_noDataValues[0];
+ else{
+ std::string errorString="Error: no valid data found";
+ throw(errorString);
+ }
}
else
return 0;
}
+//todo: what if no valid data?
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;
@@ -873,6 +1023,7 @@ template<class T> double StatFactory::cross_correlation(const std::vector<T>& x,
return sumCorrelation;
}
+//todo: nodata?
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());
@@ -884,6 +1035,7 @@ template<class T> double StatFactory::linear_regression(const std::vector<T>& x,
return (1-sumsq/var(y)/(y.size()-1));
}
+//todo: nodata?
template<class T> double StatFactory::linear_regression_err(const std::vector<T>& x, const std::vector<T>& y, double &c0, double &c1) const{
assert(x.size()==y.size());
assert(x.size());
@@ -898,6 +1050,7 @@ template<class T> double StatFactory::linear_regression_err(const std::vector<T>
//alternatively: use GNU scientific library:
// gsl_stats_correlation (const double data1[], const size_t stride1, const double data2[], const size_t stride2, const size_t n)
+//todo: nodata?
template<class T> void StatFactory::interpolateUp(const std::vector<double>& wavelengthIn, const std::vector<T>& input, const std::vector<double>& wavelengthOut, const std::string& type, std::vector<T>& output, bool verbose) const{
assert(wavelengthIn.size());
assert(input.size()==wavelengthIn.size());
@@ -969,6 +1122,7 @@ template<class T> void StatFactory::interpolateUp(const std::vector<double>& wav
// gsl_interp_accel_free(acc);
// }
+//todo: nodata?
template<class T> void StatFactory::interpolateUp(const std::vector<T>& input, std::vector<T>& output, int nbin) const
{
assert(input.size());
@@ -990,6 +1144,7 @@ template<class T> void StatFactory::interpolateUp(const std::vector<T>& input, s
}
}
+//todo: nodata?
template<class T> void StatFactory::nearUp(const std::vector<T>& input, std::vector<T>& output) const
{
assert(input.size());
@@ -1006,6 +1161,7 @@ template<class T> void StatFactory::nearUp(const std::vector<T>& input, std::vec
}
}
+//todo: nodata?
template<class T> void StatFactory::interpolateUp(double* input, int dim, std::vector<T>& output, int nbin)
{
assert(nbin);
@@ -1025,6 +1181,7 @@ template<class T> void StatFactory::interpolateUp(double* input, int dim, std::v
}
}
+//todo: nodata?
template<class T> void StatFactory::interpolateDown(const std::vector<T>& input, std::vector<T>& output, int nbin) const
{
assert(input.size());
@@ -1043,6 +1200,7 @@ template<class T> void StatFactory::interpolateDown(const std::vector<T>& input,
}
}
+//todo: nodata?
template<class T> void StatFactory::interpolateDown(double* input, int dim, std::vector<T>& output, int nbin)
{
assert(nbin);
diff --git a/src/apps/pkfilter.cc b/src/apps/pkfilter.cc
index 8c0cded..67489ad 100644
--- a/src/apps/pkfilter.cc
+++ b/src/apps/pkfilter.cc
@@ -251,6 +251,7 @@ int main(int argc,char **argv) {
for(int imask=0;imask<nodata_opt.size();++imask){
if(verbose_opt[0])
std::cout<< nodata_opt[imask] << " ";
+ filter1d.pushNoDataValue(nodata_opt[imask]);
filter2d.pushNoDataValue(nodata_opt[imask]);
}
if(verbose_opt[0])
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/pkg-grass/pktools.git
More information about the Pkg-grass-devel
mailing list