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 <kempe...@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], [kempe...@gmail.com]) +AC_INIT([pktools], [2.6.1], [kempe...@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 _______________________________________________ Pkg-grass-devel mailing list Pkg-grass-devel@lists.alioth.debian.org http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-grass-devel