[pktools] 338/375: filebuffer for ndvi in pkcomposite
Bas Couwenberg
sebastic at xs4all.nl
Wed Dec 3 21:54:28 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 629acee6e13df04d491fa2a31a68fd60ca3e2c95
Author: Pieter Kempeneers <kempenep at gmail.com>
Date: Mon Oct 13 16:27:48 2014 +0200
filebuffer for ndvi in pkcomposite
---
ChangeLog | 1 +
src/algorithms/Filter.cc | 46 ++--
src/algorithms/Filter.h | 514 +++++++++++++++++++--------------------
src/algorithms/ImgRegression.cc | 2 +-
src/apps/pkcomposite.cc | 8 +-
src/apps/pkfilter.cc | 11 +-
src/apps/pkkalman.cc | 264 ++++++++++++++++----
src/imageclasses/ImgReaderGdal.h | 27 +-
8 files changed, 514 insertions(+), 359 deletions(-)
diff --git a/ChangeLog b/ChangeLog
index 7dda69c..3eb5e35 100755
--- a/ChangeLog
+++ b/ChangeLog
@@ -310,6 +310,7 @@ version 2.5.4
Support multiple input bands when calculating statistics
- pkfilter
Support filtering and statistics in spectral domain (see ticket #43252)
+ Support different padding strategies (currently only supported for spectral/temporal filtering)
- pkextract
bug fix for proportion rule
support standard deviation rule (stdev) for polygon features (ticket #43193)
diff --git a/src/algorithms/Filter.cc b/src/algorithms/Filter.cc
index 82bd26d..4241d41 100644
--- a/src/algorithms/Filter.cc
+++ b/src/algorithms/Filter.cc
@@ -25,11 +25,13 @@ along with pktools. If not, see <http://www.gnu.org/licenses/>.
using namespace std;
filter::Filter::Filter(void)
+ : m_padding("symmetric")
{
}
filter::Filter::Filter(const vector<double> &taps)
+ : m_padding("symmetric")
{
setTaps(taps);
}
@@ -178,6 +180,7 @@ void filter::Filter::dwtCutFrom(const ImgReaderGdal& input, ImgWriterGdal& outpu
}
}
+//todo: support different padding strategies
void filter::Filter::dwtForward(std::vector<double>& data, const std::string& wavelet_type, int family){
int origsize=data.size();
//make sure data size if power of 2
@@ -196,6 +199,7 @@ void filter::Filter::dwtForward(std::vector<double>& data, const std::string& wa
gsl_wavelet_workspace_free (work);
}
+//todo: support different padding strategies
void filter::Filter::dwtInverse(std::vector<double>& data, const std::string& wavelet_type, int family){
int origsize=data.size();
//make sure data size if power of 2
@@ -213,6 +217,7 @@ void filter::Filter::dwtInverse(std::vector<double>& data, const std::string& wa
gsl_wavelet_workspace_free (work);
}
+//todo: support different padding strategies
void filter::Filter::dwtCut(std::vector<double>& data, const std::string& wavelet_type, int family, double cut){
int origsize=data.size();
//make sure data size if power of 2
@@ -241,9 +246,9 @@ void filter::Filter::dwtCut(std::vector<double>& data, const std::string& wavele
gsl_wavelet_workspace_free (work);
}
-void filter::Filter::morphology(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, int dim, short down, int offset, short verbose)
+void filter::Filter::morphology(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, int dim, short verbose)
{
- bool bverbose=(verbose>1)? true:false;
+ // bool bverbose=(verbose>1)? true:false;
Vector2d<double> lineInput(input.nrOfBand(),input.nrOfCol());
Vector2d<double> lineOutput(input.nrOfBand(),input.nrOfCol());
const char* pszMessage;
@@ -258,7 +263,8 @@ void filter::Filter::morphology(const ImgReaderGdal& input, ImgWriterGdal& outpu
vector<double> pixelOutput(input.nrOfBand());
for(int x=0;x<input.nrOfCol();++x){
pixelInput=lineInput.selectCol(x);
- morphology(pixelInput,pixelOutput,method,dim,down,offset,bverbose);
+ filter(pixelInput,pixelOutput,method,dim);
+ // morphology(pixelInput,pixelOutput,method,dim,bverbose);
for(int iband=0;iband<input.nrOfBand();++iband)
lineOutput[iband][x]=pixelOutput[iband];
}
@@ -275,13 +281,13 @@ void filter::Filter::morphology(const ImgReaderGdal& input, ImgWriterGdal& outpu
}
}
-void filter::Filter::smooth(const ImgReaderGdal& input, ImgWriterGdal& output, short dim, short down, int offset)
+void filter::Filter::smooth(const ImgReaderGdal& input, ImgWriterGdal& output, short dim)
{
assert(dim>0);
m_taps.resize(dim);
for(int itap=0;itap<dim;++itap)
m_taps[itap]=1.0/dim;
- filter(input,output,down,offset);
+ filter(input,output);
}
// void filter::Filter::smoothnodata(const ImgReaderGdal& input, ImgWriterGdal& output, short dim, short down, int offset)
@@ -293,7 +299,7 @@ void filter::Filter::smooth(const ImgReaderGdal& input, ImgWriterGdal& output, s
// filter(input,output,down,offset);
// }
-void filter::Filter::filter(const ImgReaderGdal& input, ImgWriterGdal& output, short down, int offset)
+void filter::Filter::filter(const ImgReaderGdal& input, ImgWriterGdal& output)
{
Vector2d<double> lineInput(input.nrOfBand(),input.nrOfCol());
Vector2d<double> lineOutput(input.nrOfBand(),input.nrOfCol());
@@ -309,7 +315,7 @@ void filter::Filter::filter(const ImgReaderGdal& input, ImgWriterGdal& output, s
vector<double> pixelOutput(input.nrOfBand());
for(int x=0;x<input.nrOfCol();++x){
pixelInput=lineInput.selectCol(x);
- filter(pixelInput,pixelOutput,down,offset);
+ filter(pixelInput,pixelOutput);
for(int iband=0;iband<input.nrOfBand();++iband)
lineOutput[iband][x]=pixelOutput[iband];
}
@@ -326,10 +332,10 @@ void filter::Filter::filter(const ImgReaderGdal& input, ImgWriterGdal& output, s
}
}
-void filter::Filter::stat(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, short down, int offset)
+void filter::Filter::stat(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method)
{
Vector2d<double> lineInput(input.nrOfBand(),input.nrOfCol());
- assert(output.nrOfCol()==(input.nrOfCol()+down-1)/down);
+ assert(output.nrOfCol()==input.nrOfCol());
vector<double> lineOutput(output.nrOfCol());
statfactory::StatFactory stat;
const char* pszMessage;
@@ -338,36 +344,32 @@ void filter::Filter::stat(const ImgReaderGdal& input, ImgWriterGdal& output, con
double progress=0;
pfnProgress(progress,pszMessage,pProgressArg);
for(int y=0;y<input.nrOfRow();++y){
- if((y+1+down/2)%down)
- continue;
for(int iband=0;iband<input.nrOfBand();++iband)
input.readData(lineInput[iband],GDT_Float64,y,iband);
vector<double> pixelInput(input.nrOfBand());
for(int x=0;x<input.nrOfCol();++x){
- if((x+1+down/2)%down)
- continue;
pixelInput=lineInput.selectCol(x);
switch(getFilterType(method)){
case(filter::median):
- lineOutput[(x-offset+down-1)/down]=stat.median(pixelInput);
+ lineOutput[x]=stat.median(pixelInput);
break;
case(filter::min):
- lineOutput[(x-offset+down-1)/down]=stat.mymin(pixelInput);
+ lineOutput[x]=stat.mymin(pixelInput);
break;
case(filter::max):
- lineOutput[(x-offset+down-1)/down]=stat.mymax(pixelInput);
+ lineOutput[x]=stat.mymax(pixelInput);
break;
case(filter::sum):
- lineOutput[(x-offset+down-1)/down]=stat.sum(pixelInput);
+ lineOutput[x]=stat.sum(pixelInput);
break;
case(filter::var):
- lineOutput[(x-offset+down-1)/down]=stat.var(pixelInput);
+ lineOutput[x]=stat.var(pixelInput);
break;
case(filter::stdev):
- lineOutput[(x-offset+down-1)/down]=sqrt(stat.var(pixelInput));
+ lineOutput[x]=sqrt(stat.var(pixelInput));
break;
case(filter::mean):
- lineOutput[(x-offset+down-1)/down]=stat.mean(pixelInput);
+ lineOutput[x]=stat.mean(pixelInput);
break;
default:
std::string errorString="method not supported";
@@ -376,10 +378,10 @@ void filter::Filter::stat(const ImgReaderGdal& input, ImgWriterGdal& output, con
}
}
try{
- output.writeData(lineOutput,GDT_Float64,y/down);
+ output.writeData(lineOutput,GDT_Float64,y);
}
catch(string errorstring){
- cerr << errorstring << "in line " << y/down << endl;
+ cerr << errorstring << "in line " << y << endl;
}
progress=(1.0+y)/output.nrOfRow();
pfnProgress(progress,pszMessage,pProgressArg);
diff --git a/src/algorithms/Filter.h b/src/algorithms/Filter.h
index 70955a1..9b12410 100644
--- a/src/algorithms/Filter.h
+++ b/src/algorithms/Filter.h
@@ -35,12 +35,19 @@ namespace filter
enum FILTER_TYPE { median=0, var=1 , min=2, max=3, sum=4, mean=5, minmax=6, dilate=7, erode=8, close=9, open=10, homog=11, sobelx=12, sobely=13, sobelxy=14, sobelyx=-14, smooth=15, density=16, majority=17, mixed=18, smoothnodata=19, threshold=20, ismin=21, ismax=22, heterog=23, order=24, stdev=25, dwt=26, dwti=27, dwt_cut=28, dwt_cut_from=29};
+ enum PADDING { symmetric=0, replicate=1, circular=2, constant=3};
+
class Filter
{
public:
Filter(void);
Filter(const std::vector<double> &taps);
virtual ~Filter(){};
+
+ void setPadding(const std::string& padString){
+ m_padding=padString;
+ };
+
static const gsl_wavelet_type* getWaveletType(const std::string waveletType){
if(waveletType=="daubechies") return(gsl_wavelet_daubechies);
if(waveletType=="daubechies_centered") return(gsl_wavelet_daubechies_centered);
@@ -54,20 +61,21 @@ public:
initFilterMap(m_filterMap);
return m_filterMap[filterType];
};
+
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);};
- template<class T> void filter(const std::vector<T>& input, std::vector<T>& output, int down=1, int offset=0);
+ 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, int down=1, int offset=0);
- template<class T> void filter(T* input, int inputSize, std::vector<T>& output, int down=1, int offset=0);
- template<class T> void smooth(T* input, int inputSize, std::vector<T>& output, short dim, int down=1, int offset=0);
- 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=false);
- void morphology(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, int dim, short down=1, int offset=0, short verbose=0);
- void filter(const ImgReaderGdal& input, ImgWriterGdal& output, short down=1, int offset=0);
- void stat(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, short down=1, int offset=0);
+ template<class T> void smooth(const std::vector<T>& input, std::vector<T>& output, short dim);
+ template<class T> void filter(T* input, int inputSize, std::vector<T>& output);
+ template<class T> void smooth(T* input, int inputSize, std::vector<T>& output, short dim);
+ //template<class T> void morphology(const std::vector<T>& input, std::vector<T>& output, const std::string& method, int dim, bool verbose=false);
+ void morphology(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, int dim, short verbose=0);
+ void filter(const ImgReaderGdal& input, ImgWriterGdal& output);
+ void stat(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method);
void filter(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, int dim);
- void smooth(const ImgReaderGdal& input, ImgWriterGdal& output, short dim, short down=1, int offset=0);
+ void smooth(const ImgReaderGdal& input, ImgWriterGdal& output, short dim);
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);
@@ -118,11 +126,24 @@ private:
m_filterMap["order"]=filter::order;
m_filterMap["median"]=filter::median;
}
+
+
+ static PADDING getPadding(const std::string& padString){
+ std::map<std::string, PADDING> padMap;
+ padMap["constant"]=filter::constant;
+ padMap["symmetric"]=filter::symmetric;
+ padMap["replicate"]=filter::replicate;
+ padMap["circular"]=filter::circular;
+ return(padMap[padString]);
+ };
+
std::vector<double> m_taps;
std::vector<short> m_class;
std::vector<short> m_mask;
+ std::string m_padding;
};
+
//input[band], output
//returns wavelength for which srf is maximum
template<class T> double Filter::applySrf(const std::vector<double> &wavelengthIn, const std::vector<T>& input, const Vector2d<double>& srf, const std::string& interpolationType, T& output, double delta, bool normalize, bool verbose)
@@ -195,15 +216,6 @@ private:
gsl_spline_free(splineOut);
gsl_interp_accel_free(accOut);
- // double maxResponse=0;
- // int maxIndex=0;
- // for(int index=0;index<srf[1].size();++index){
- // if(maxResponse<srf[1][index]){
- // maxResponse=srf[1][index];
- // maxIndex=index;
- // }
- // }
- // return(srf[0][maxIndex]);
return(centreWavelength);
}
@@ -291,15 +303,6 @@ private:
gsl_spline_free(splineOut);
gsl_interp_accel_free(accOut);
- // double maxResponse=0;
- // int maxIndex=0;
- // for(int index=0;index<srf[1].size();++index){
- // if(maxResponse<srf[1][index]){
- // maxResponse=srf[1][index];
- // maxIndex=index;
- // }
- // }
- // return(srf[0][maxIndex]);
return(centreWavelength);
}
@@ -405,51 +408,88 @@ template<class T> void Filter::applyFwhm(const std::vector<double> &wavelengthIn
}
}
- template<class T> void Filter::smooth(const std::vector<T>& input, std::vector<T>& output, short dim, int down, int offset)
+ template<class T> void Filter::smooth(const std::vector<T>& input, std::vector<T>& output, short dim)
{
assert(dim>0);
m_taps.resize(dim);
for(int itap=0;itap<dim;++itap)
m_taps[itap]=1.0/dim;
- filter(input,output,down,offset);
+ filter(input,output);
}
-template<class T> void Filter::filter(const std::vector<T>& input, std::vector<T>& output, int down, int offset)
+template<class T> void Filter::filter(const std::vector<T>& input, std::vector<T>& output)
{
- output.resize((input.size()-offset+down-1)/down);
+ assert(input.size()>m_taps.size());
+ output.resize(input.size());
int i=0;
- //start: extend input with mirrored version of itself
- for(i=offset;i<m_taps.size()/2;++i){
- if((i-offset)%down)
- continue;
+ //start: extend input by padding
+ for(i=0;i<m_taps.size()/2;++i){
//todo:introduce nodata
- output[(i-offset+down-1)/down]=m_taps[m_taps.size()/2]*input[i];
- for(int t=1;t<=m_taps.size()/2;++t)
- output[(i-offset+down-1)/down]+=(m_taps[m_taps.size()/2+t]+m_taps[m_taps.size()/2-t])*input[i+t];
+ 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];
+ if(i>=t)
+ output[i]+=m_taps[m_taps.size()/2-t]*input[i-t];
+ else{
+ switch(getPadding(m_padding)){
+ case(replicate):
+ output[i]+=m_taps[m_taps.size()/2-t]*input[0];
+ break;
+ case(circular):
+ output[i]+=m_taps[m_taps.size()/2-t]*input[input.size()+i-t];
+ break;
+ case(constant):
+ output[i]+=m_taps[m_taps.size()/2-t]*0;
+ break;
+ case(symmetric):
+ default:
+ output[i]+=m_taps[m_taps.size()/2-t]*input[t-i];
+ break;
+ }
+ }
+ }
}
//main
- for(i=offset+m_taps.size()/2;i<input.size()-m_taps.size()/2;++i){
- if((i-offset)%down)
- continue;
+ for(i=m_taps.size()/2;i<input.size()-m_taps.size()/2;++i){
//todo:introduce nodata
T leaveOut=(*(m_taps.begin()))*input[i-m_taps.size()/2];
T include=(m_taps.back())*input[i+m_taps.size()/2];
- output[(i-offset+down-1)/down]=0;
+ output[i]=0;
for(int t=0;t<m_taps.size();++t)
- output[(i-offset+down-1)/down]+=input[i-m_taps.size()/2+t]*m_taps[t];
+ output[i]+=input[i-m_taps.size()/2+t]*m_taps[t];
}
- //end: extend input with mirrored version of itself
+ //end: extend input by padding
for(i=input.size()-m_taps.size()/2;i<input.size();++i){
- if((i-offset)%down)
- continue;
//todo:introduce nodata
- output[(i-offset+down-1)/down]=m_taps[m_taps.size()/2]*input[i];
+ output[i]=m_taps[m_taps.size()/2]*input[i];
//todo:introduce nodata
- for(int t=1;t<=m_taps.size()/2;++t)
- output[(i-offset+down-1)/down]+=(m_taps[m_taps.size()/2+t]+m_taps[m_taps.size()/2-t])*input[i-t];
+ 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())
+ output[i]+=m_taps[m_taps.size()/2+t]*input[i+t];
+ else{
+ switch(getPadding(m_padding)){
+ case(replicate):
+ output[i]+=m_taps[m_taps.size()/2+t]*input.back();
+ break;
+ case(circular):
+ output[i]+=m_taps[m_taps.size()/2+t]*input[t-1];
+ break;
+ case(constant):
+ output[i]+=m_taps[m_taps.size()/2+t]*0;
+ break;
+ case(symmetric):
+ default:
+ output[i]+=m_taps[m_taps.size()/2+t]*input[i-t];
+ break;
+ }
+ }
+ //output[i]+=(m_taps[m_taps.size()/2+t]+m_taps[m_taps.size()/2-t])*input[i-t];
+ }
}
}
+//todo: filling statBuffer can be optimized (no need to clear and fill entire buffer, just push back new value...)
template<class T> void Filter::filter(const std::vector<T>& input, std::vector<T>& output, const std::string& method, int dim)
{
bool verbose=false;
@@ -459,7 +499,7 @@ template<class T> void Filter::filter(const std::vector<T>& input, std::vector<T
statfactory::StatFactory stat;
std::vector<T> statBuffer;
short binValue=0;
- //start: extend input with mirrored version of itself
+ //start: extend input by padding
for(i=0;i<dim/2;++i){
binValue=0;
for(int iclass=0;iclass<m_class.size();++iclass){
@@ -472,36 +512,62 @@ template<class T> void Filter::filter(const std::vector<T>& input, std::vector<T
statBuffer.push_back(binValue);
else
statBuffer.push_back(input[i]);
+
for(int t=1;t<=dim/2;++t){
- binValue=0;
+ T theValue=input[i+t];
for(int iclass=0;iclass<m_class.size();++iclass){
- if(input[i+t]==m_class[iclass]){
+ if(theValue==m_class[iclass]){
binValue=m_class[0];
break;
}
}
- if(m_class.size()){
- statBuffer.push_back(binValue);
- statBuffer.push_back(binValue);
+ if(m_class.size())
+ statBuffer.push_back(binValue);
+ else
+ statBuffer.push_back(theValue);
+
+ if(i>=t){
+ theValue=input[i-t];
}
else{
- statBuffer.push_back(input[i+t]);
- statBuffer.push_back(input[i+t]);
+ switch(getPadding(m_padding)){
+ case(replicate):
+ theValue=input[0];
+ break;
+ case(circular):
+ theValue=input[input.size()+i-t];
+ break;
+ case(constant):
+ theValue=0;
+ break;
+ case(symmetric):
+ default:
+ theValue=input[t-i];
+ break;
+ }
+ }
+ for(int iclass=0;iclass<m_class.size();++iclass){
+ if(theValue==m_class[iclass]){
+ binValue=m_class[0];
+ break;
+ }
}
+ if(m_class.size())
+ statBuffer.push_back(binValue);
+ else
+ statBuffer.push_back(theValue);
}
- /* assert(statBuffer.size()==dim); */
- /* if((i-offset)%down){ */
- /* statBuffer.clear(); */
- /* continue; */
- /* } */
+
switch(getFilterType(method)){
case(filter::median):
output[i]=stat.median(statBuffer);
break;
case(filter::min):
+ case(filter::erode):
output[i]=stat.mymin(statBuffer);
break;
case(filter::max):
+ case(filter::dilate):
output[i]=stat.mymax(statBuffer);
break;
case(filter::sum):
@@ -535,15 +601,16 @@ template<class T> void Filter::filter(const std::vector<T>& input, std::vector<T
else
statBuffer.push_back(input[i-dim/2+t]);
}
- /* assert(statBuffer.size()==dim); */
switch(getFilterType(method)){
case(filter::median):
output[i]=stat.median(statBuffer);
break;
case(filter::min):
+ case(filter::erode):
output[i]=stat.mymin(statBuffer);
break;
case(filter::max):
+ case(filter::dilate):
output[i]=stat.mymax(statBuffer);
break;
case(filter::sum):
@@ -562,262 +629,173 @@ template<class T> void Filter::filter(const std::vector<T>& input, std::vector<T
}
statBuffer.clear();
}
- //end: extend input with mirrored version of itself
+ //end: extend input by padding
for(i=input.size()-dim/2;i<input.size();++i){
- binValue=0;
- for(int iclass=0;iclass<m_class.size();++iclass){
- if(input[i]==m_class[iclass]){
- binValue=m_class[0];
- break;
- }
- }
- if(m_class.size())
- statBuffer.push_back(binValue);
- else
- statBuffer.push_back(input[i]);
- for(int t=1;t<=dim/2;++t){
- binValue=0;
- for(int iclass=0;iclass<m_class.size();++iclass){
- if(input[i-t]==m_class[iclass]){
- binValue=m_class[0];
- break;
- }
- }
- if(m_class.size()){
- statBuffer.push_back(binValue);
- statBuffer.push_back(binValue);
- }
- else{
- statBuffer.push_back(input[i-t]);
- statBuffer.push_back(input[i-t]);
- }
- }
- switch(getFilterType(method)){
- case(filter::median):
- output[i]=stat.median(statBuffer);
- break;
- case(filter::min):
- output[i]=stat.mymin(statBuffer);
- break;
- case(filter::max):
- output[i]=stat.mymax(statBuffer);
- break;
- case(filter::sum):
- output[i]=sqrt(stat.sum(statBuffer));
- break;
- case(filter::var):
- output[i]=stat.var(statBuffer);
- break;
- case(filter::mean):
- output[i]=stat.mean(statBuffer);
- break;
- default:
- std::string errorString="method not supported";
- throw(errorString);
- break;
- }
- }
-}
-
- //todo: this function is redundant, can be incorporated within filter
-template<class T> void Filter::morphology(const std::vector<T>& input, std::vector<T>& output, const std::string& method, int dim, short down, int offset, bool verbose)
-{
- assert(dim);
- output.resize((input.size()-offset+down-1)/down);
- int i=0;
- statfactory::StatFactory stat;
- std::vector<T> statBuffer;
- short binValue=0;
- //start: extend input with mirrored version of itself
- for(i=offset;i<dim/2;++i){
binValue=0;
for(int iclass=0;iclass<m_class.size();++iclass){
if(input[i]==m_class[iclass]){
- binValue=m_class[0];
- break;
+ binValue=m_class[0];
+ break;
}
}
if(m_class.size())
statBuffer.push_back(binValue);
else
statBuffer.push_back(input[i]);
+
for(int t=1;t<=dim/2;++t){
- binValue=0;
+ T theValue=input[i-t];
for(int iclass=0;iclass<m_class.size();++iclass){
- if(input[i+t]==m_class[iclass]){
- binValue=m_class[0];
- break;
- }
- }
- if(m_class.size()){
- statBuffer.push_back(binValue);
- statBuffer.push_back(binValue);
+ if(theValue==m_class[iclass]){
+ binValue=m_class[0];
+ break;
+ }
}
+ if(m_class.size())
+ statBuffer.push_back(binValue);
+ else
+ statBuffer.push_back(theValue);
+ if(i+t<input.size())
+ theValue=input[i+t];
else{
- statBuffer.push_back(input[i+t]);
- statBuffer.push_back(input[i+t]);
+ switch(getPadding(m_padding)){
+ case(replicate):
+ theValue=input.back();
+ break;
+ case(circular):
+ theValue=input[t-1];
+ break;
+ case(constant):
+ theValue=0;
+ break;
+ case(symmetric):
+ default:
+ theValue=input[i-t];
+ break;
+ }
}
- }
- /* assert(statBuffer.size()==dim); */
- if((i-offset)%down){
- statBuffer.clear();
- continue;
- }
- switch(getFilterType(method)){
- case(filter::dilate):
- output[(i-offset+down-1)/down]=stat.mymax(statBuffer);
- break;
- case(filter::erode):
- output[(i-offset+down-1)/down]=stat.mymin(statBuffer);
- break;
- default:
- std::string errorString="method not supported";
- throw(errorString);
- break;
- }
- if(verbose){
- std::cout << "buffer: ";
- for(int ibuf=0;ibuf<statBuffer.size();++ibuf)
- std::cout << statBuffer[ibuf] << " ";
- std::cout << "->" << output[(i-offset+down-1)/down] << std::endl;
- }
- }
- //main
- statBuffer.clear();
- for(i=offset+dim/2;i<input.size()-dim/2;++i){
- binValue=0;
- for(int t=0;t<dim;++t){
for(int iclass=0;iclass<m_class.size();++iclass){
- if(input[i-dim/2+t]==m_class[iclass]){
- binValue=m_class[0];
- break;
- }
+ if(theValue==m_class[iclass]){
+ binValue=m_class[0];
+ break;
+ }
}
if(m_class.size())
- statBuffer.push_back(binValue);
+ statBuffer.push_back(binValue);
else
- statBuffer.push_back(input[i-dim/2+t]);
- }
- /* assert(statBuffer.size()==dim); */
- if((i-offset)%down){
- statBuffer.clear();
- continue;
+ statBuffer.push_back(theValue);
}
switch(getFilterType(method)){
- case(filter::dilate):
- output[(i-offset+down-1)/down]=stat.mymax(statBuffer);
+ case(filter::median):
+ output[i]=stat.median(statBuffer);
break;
+ case(filter::min):
case(filter::erode):
- output[(i-offset+down-1)/down]=stat.mymin(statBuffer);
- break;
- default:
- std::string errorString="method not supported";
- throw(errorString);
+ output[i]=stat.mymin(statBuffer);
break;
- }
- if(verbose){
- std::cout << "buffer: ";
- for(int ibuf=0;ibuf<statBuffer.size();++ibuf)
- std::cout << statBuffer[ibuf] << " ";
- std::cout << "->" << output[(i-offset+down-1)/down] << std::endl;
- }
- statBuffer.clear();
- }
- //end: extend input with mirrored version of itself
- for(i=input.size()-dim/2;i<input.size();++i){
- binValue=0;
- for(int iclass=0;iclass<m_class.size();++iclass){
- if(input[i]==m_class[iclass]){
- binValue=m_class[0];
- break;
- }
- }
- if(m_class.size())
- statBuffer.push_back(binValue);
- else
- statBuffer.push_back(input[i]);
- for(int t=1;t<=dim/2;++t){
- binValue=0;
- for(int iclass=0;iclass<m_class.size();++iclass){
- if(input[i-t]==m_class[iclass]){
- binValue=m_class[0];
- break;
- }
- }
- if(m_class.size()){
- statBuffer.push_back(binValue);
- statBuffer.push_back(binValue);
- }
- else{
- statBuffer.push_back(input[i-t]);
- statBuffer.push_back(input[i-t]);
- }
- }
- if((i-offset)%down){
- statBuffer.clear();
- continue;
- }
- switch(getFilterType(method)){
+ case(filter::max):
case(filter::dilate):
- output[(i-offset+down-1)/down]=stat.mymax(statBuffer);
+ output[i]=stat.mymax(statBuffer);
break;
- case(filter::erode):
- output[(i-offset+down-1)/down]=stat.mymin(statBuffer);
+ case(filter::sum):
+ output[i]=sqrt(stat.sum(statBuffer));
+ break;
+ case(filter::var):
+ output[i]=stat.var(statBuffer);
+ break;
+ case(filter::mean):
+ output[i]=stat.mean(statBuffer);
break;
default:
std::string errorString="method not supported";
throw(errorString);
break;
}
- if(verbose){
- std::cout << "buffer: ";
- for(int ibuf=0;ibuf<statBuffer.size();++ibuf)
- std::cout << statBuffer[ibuf] << " ";
- std::cout << "->" << output[(i-offset+down-1)/down] << std::endl;
- }
}
-}
+ }
- template<class T> void Filter::smooth(T* input, int inputSize, std::vector<T>& output, short dim, int down, int offset)
+ template<class T> void Filter::smooth(T* input, int inputSize, std::vector<T>& output, short dim)
{
assert(dim>0);
m_taps.resize(dim);
for(int itap=0;itap<dim;++itap)
m_taps[itap]=1.0/dim;
- filter(input,output,down,offset);
+ filter(input,output);
}
-template<class T> void Filter::filter(T* input, int inputSize, std::vector<T>& output, int down, int offset)
+template<class T> void Filter::filter(T* input, int inputSize, std::vector<T>& output)
{
- output.resize((inputSize-offset+down-1)/down);
+ assert(inputSize>m_taps.size());
+ output.resize(inputSize);
int i=0;
- //start: extend input with mirrored version of itself
- for(i=offset;i<m_taps.size()/2;++i){
- if((i-offset)%down)
- continue;
- output[(i-offset+down-1)/down]=m_taps[m_taps.size()/2]*input[i];
- for(int t=1;t<=m_taps.size()/2;++t)
- output[(i-offset+down-1)/down]+=(m_taps[m_taps.size()/2+t]+m_taps[m_taps.size()/2-t])*input[i+t];
+
+ //start: extend input by padding
+ for(i=0;i<m_taps.size()/2;++i){
+ //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];
+ if(i>=t)
+ output[i]+=m_taps[m_taps.size()/2-t]*input[i-t];
+ else{
+ switch(getPadding(m_padding)){
+ case(replicate):
+ output[i]+=m_taps[m_taps.size()/2-t]*input[0];
+ break;
+ case(circular):
+ output[i]+=m_taps[m_taps.size()/2-t]*input[input.size()+i-t];
+ break;
+ case(constant):
+ output[i]+=m_taps[m_taps.size()/2-t]*0;
+ break;
+ case(symmetric):
+ default:
+ output[i]+=m_taps[m_taps.size()/2-t]*input[t-i];
+ break;
+ }
+ }
+ }
}
//main
- for(i=offset+m_taps.size()/2;i<inputSize-m_taps.size()/2;++i){
- if((i-offset)%down)
- continue;
+ for(i=m_taps.size()/2;i<input.size()-m_taps.size()/2;++i){
+ //todo:introduce nodata
T leaveOut=(*(m_taps.begin()))*input[i-m_taps.size()/2];
T include=(m_taps.back())*input[i+m_taps.size()/2];
- output[(i-offset+down-1)/down]=0;
+ output[i]=0;
for(int t=0;t<m_taps.size();++t)
- output[(i-offset+down-1)/down]+=input[i-m_taps.size()/2+t]*m_taps[t];
+ output[i]+=input[i-m_taps.size()/2+t]*m_taps[t];
}
- //end: extend input with mirrored version of itself
- for(i=inputSize-m_taps.size()/2;i<inputSize;++i){
- if((i-offset)%down)
- continue;
- output[(i-offset+down-1)/down]=m_taps[m_taps.size()/2]*input[i];
- for(int t=1;t<=m_taps.size()/2;++t)
- output[(i-offset+down-1)/down]+=(m_taps[m_taps.size()/2+t]+m_taps[m_taps.size()/2-t])*input[i-t];
+ //end: extend input by padding
+ for(i=input.size()-m_taps.size()/2;i<input.size();++i){
+ //todo:introduce nodata
+ output[i]=m_taps[m_taps.size()/2]*input[i];
+ //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())
+ output[i]+=m_taps[m_taps.size()/2+t]*input[i+t];
+ else{
+ switch(getPadding(m_padding)){
+ case(replicate):
+ output[i]+=m_taps[m_taps.size()/2+t]*input.back();
+ break;
+ case(circular):
+ output[i]+=m_taps[m_taps.size()/2+t]*input[t-1];
+ break;
+ case(constant):
+ output[i]+=m_taps[m_taps.size()/2+t]*0;
+ break;
+ case(symmetric):
+ default:
+ output[i]+=m_taps[m_taps.size()/2+t]*input[i-t];
+ break;
+ }
+ }
+ }
}
}
+
}
#endif /* _MYFILTER_H_ */
diff --git a/src/algorithms/ImgRegression.cc b/src/algorithms/ImgRegression.cc
index 55b3c2f..e1de847 100644
--- a/src/algorithms/ImgRegression.cc
+++ b/src/algorithms/ImgRegression.cc
@@ -81,7 +81,7 @@ double ImgRegression::getRMSE(const ImgReaderGdal& imgReader1, const ImgReaderGd
}
}
double err=0;
- if(buffer1.size()||buffer2.size()){
+ if(buffer1.size()&&buffer2.size()){
statfactory::StatFactory stat;
err=stat.linear_regression_err(buffer1,buffer2,c0,c1);
}
diff --git a/src/apps/pkcomposite.cc b/src/apps/pkcomposite.cc
index f53135a..0854009 100644
--- a/src/apps/pkcomposite.cc
+++ b/src/apps/pkcomposite.cc
@@ -627,8 +627,8 @@ int main(int argc, char *argv[])
val_new=(readCol-0.5-lowerCol)*readBuffer[iband][upperCol-startCol]+(1-readCol+0.5+lowerCol)*readBuffer[iband][lowerCol-startCol];
writeBuffer[iband][ib]=val_new;
}
- // fileBuffer[ib]=ifile;
- ++fileBuffer[ib];
+ fileBuffer[ib]=ifile;
+ // ++fileBuffer[ib];
}
break;
default:
@@ -642,8 +642,8 @@ int main(int argc, char *argv[])
val_new=readBuffer[iband][readCol-startCol];
writeBuffer[iband][ib]=val_new;
}
- ++fileBuffer[ib];
- // fileBuffer[ib]=ifile;
+ fileBuffer[ib]=ifile;
+ // ++fileBuffer[ib];
}
break;
}
diff --git a/src/apps/pkfilter.cc b/src/apps/pkfilter.cc
index 4018896..21d43cc 100644
--- a/src/apps/pkfilter.cc
+++ b/src/apps/pkfilter.cc
@@ -54,6 +54,7 @@ int main(int argc,char **argv) {
Optionpk<short> nodata_opt("nodata", "nodata", "nodata value(s) for smoothnodata filter");
Optionpk<std::string> tap_opt("tap", "tap", "text file containing taps used for spatial filtering (from ul to lr). Use dimX and dimY to specify tap dimensions in x and y. Leave empty for not using taps");
Optionpk<double> tapz_opt("tapz", "tapz", "taps used for spectral filtering");
+ Optionpk<string> padding_opt("pad","pad", "Padding method for filtering (how to handle edge effects). Choose between: symmetric, replicate, circular, constant (pad with 0).", "symmetric");
Optionpk<double> fwhm_opt("fwhm", "fwhm", "list of full width half to apply spectral filtering (-fwhm band1 -fwhm band2 ...)");
Optionpk<std::string> srf_opt("srf", "srf", "list of ASCII files containing spectral response functions (two columns: wavelength response)");
Optionpk<double> wavelengthIn_opt("win", "wavelengthIn", "list of wavelengths in input spectrum (-win band1 -win band2 ...)");
@@ -92,6 +93,7 @@ int main(int argc,char **argv) {
nodata_opt.retrieveOption(argc,argv);
tap_opt.retrieveOption(argc,argv);
tapz_opt.retrieveOption(argc,argv);
+ padding_opt.retrieveOption(argc,argv);
fwhm_opt.retrieveOption(argc,argv);
srf_opt.retrieveOption(argc,argv);
wavelengthIn_opt.retrieveOption(argc,argv);
@@ -222,6 +224,9 @@ int main(int argc,char **argv) {
filter2d::Filter2d filter2d;
filter::Filter filter1d;
+ if(verbose_opt[0])
+ cout << "Set padding to " << padding_opt[0] << endl;
+ filter1d.setPadding(padding_opt[0]);
if(class_opt.size()){
if(verbose_opt[0])
std::cout<< "class values: ";
@@ -278,7 +283,7 @@ int main(int argc,char **argv) {
std::cout<< std::endl;
}
filter1d.setTaps(tapz_opt);
- filter1d.filter(input,output,down_opt[0]);
+ filter1d.filter(input,output);
}
else if(fwhm_opt.size()){
if(verbose_opt[0])
@@ -385,7 +390,7 @@ int main(int argc,char **argv) {
if(dimZ_opt.size()){
if(verbose_opt[0])
std::cout<< "1-D filtering: dilate" << std::endl;
- filter1d.morphology(input,output,"dilate",dimZ_opt[0],1,0,verbose_opt[0]);
+ filter1d.morphology(input,output,"dilate",dimZ_opt[0],verbose_opt[0]);
}
else{
filter2d.morphology(input,output,"dilate",dimX_opt[0],dimY_opt[0],angle_opt,disc_opt[0]);
@@ -710,7 +715,7 @@ int main(int argc,char **argv) {
default:
if(dimZ_opt.size()){
if(dimZ_opt[0]==1)
- filter1d.stat(input,output,method_opt[0],down_opt[0]);
+ filter1d.stat(input,output,method_opt[0]);
else{
assert(down_opt[0]==1);//not implemented yet...
filter1d.filter(input,output,method_opt[0],dimZ_opt[0]);
diff --git a/src/apps/pkkalman.cc b/src/apps/pkkalman.cc
index eff5c4e..04daa31 100644
--- a/src/apps/pkkalman.cc
+++ b/src/apps/pkkalman.cc
@@ -41,7 +41,6 @@ int main(int argc,char **argv) {
Optionpk<string> outputfw_opt("ofw", "outputfw", "Output raster dataset for forward model");
Optionpk<string> outputbw_opt("obw", "outputbw", "Output raster dataset for backward model");
Optionpk<string> outputfb_opt("ofb", "outputfb", "Output raster dataset for smooth model");
- Optionpk<float> threshold_opt("th", "threshold", "threshold for selecting samples (randomly). Provide probability in percentage (>0) or absolute (<0).", 0);
Optionpk<double> modnodata_opt("modnodata", "modnodata", "invalid value for model input", 0);
Optionpk<double> obsnodata_opt("obsnodata", "obsnodata", "invalid value for observation input", 0);
Optionpk<double> modoffset_opt("modoffset", "modoffset", "offset used to read model input dataset (value=offset+scale*readValue", 0);
@@ -55,6 +54,9 @@ int main(int argc,char **argv) {
// Optionpk<double> regTime_opt("rt", "regtime", "Relative Weight for regression in time series", 1.0);
// Optionpk<double> regSensor_opt("rs", "regsensor", "Relative Weight for regression in sensor series", 1.0);
Optionpk<int> down_opt("down", "down", "Downsampling factor for reading model data to calculate regression", 9);
+ Optionpk<float> threshold_opt("th", "threshold", "threshold for selecting samples (randomly). Provide probability in percentage (>0) or absolute (<0).", 0);
+ Optionpk<int> minreg_opt("minreg", "minreg", "Minimum number of pixels to take into account for regression", 5, 2);
+ Optionpk<unsigned short> window_opt("win", "window", "window size for calculating regression (use 0 for global)", 0);
Optionpk<string> oformat_opt("of", "oformat", "Output image format (see also gdal_translate). Empty string: inherit from input image","GTiff",2);
Optionpk<string> option_opt("co", "co", "Creation option for output file. Multiple options can be specified.");
Optionpk<short> verbose_opt("v", "verbose", "verbose mode when positive", 0);
@@ -70,7 +72,6 @@ int main(int argc,char **argv) {
outputfw_opt.retrieveOption(argc,argv);
outputbw_opt.retrieveOption(argc,argv);
outputfb_opt.retrieveOption(argc,argv);
- threshold_opt.retrieveOption(argc,argv);
modnodata_opt.retrieveOption(argc,argv);
obsnodata_opt.retrieveOption(argc,argv);
modoffset_opt.retrieveOption(argc,argv);
@@ -84,6 +85,9 @@ int main(int argc,char **argv) {
// regTime_opt.retrieveOption(argc,argv);
// regSensor_opt.retrieveOption(argc,argv);
down_opt.retrieveOption(argc,argv);
+ threshold_opt.retrieveOption(argc,argv);
+ minreg_opt.retrieveOption(argc,argv);
+ window_opt.retrieveOption(argc,argv);
oformat_opt.retrieveOption(argc,argv);
option_opt.retrieveOption(argc,argv);
verbose_opt.retrieveOption(argc,argv);
@@ -156,6 +160,7 @@ int main(int argc,char **argv) {
}
statfactory::StatFactory stat;
+ stat.setNoDataValues(modnodata_opt);
imgregression::ImgRegression imgreg;
// vector<ImgReaderGdal> imgReaderModel(model_opt.size());
// vector<ImgReaderGdal> imgReaderObs(observation_opt.size());
@@ -194,6 +199,11 @@ int main(int argc,char **argv) {
imgreg.setDown(down_opt[0]);
imgreg.setThreshold(threshold_opt[0]);
+ double c0modGlobal=0;//time regression coefficient c0 calculated on entire image
+ double c1modGlobal=0;//time regression coefficient c1 calculated on entire image
+ double c0mod=0;//time regression coefficient c0 calculated on local window
+ double c1mod=0;//time regression coefficient c1 calculated on local window
+
double c0obs=0;
double c1obs=0;
double errObs=uncertNodata_opt[0];//start with high initial value in case we do not have first observation at time 0
@@ -403,14 +413,12 @@ int main(int argc,char **argv) {
//to keep it general, we must redo it (overlap might have changed)
pfnProgress(progress,pszMessage,pProgressArg);
- double c0mod=0;
- double c1mod=0;
if(verbose_opt[0])
cout << "Calculating regression for " << imgReaderModel1.getFileName() << " " << imgReaderModel2.getFileName() << endl;
- double errMod=imgreg.getRMSE(imgReaderModel1,imgReaderModel2,c0mod,c1mod);
- // double errMod=imgreg.getRMSE(imgReaderModel1,imgReaderModel2,c0mod,c1mod,verbose_opt[0]);
+ double errMod=imgreg.getRMSE(imgReaderModel1,imgReaderModel2,c0modGlobal,c1modGlobal);
+ // double errMod=imgreg.getRMSE(imgReaderModel1,imgReaderModel2,c0modGlobal,c1modGlobal,verbose_opt[0]);
bool update=false;
if(obsindex<relobsindex.size()){
@@ -446,50 +454,125 @@ int main(int argc,char **argv) {
imgReaderEst.setOffset(obsoffset_opt[0]);
imgReaderEst.setScale(obsscale_opt[0]);
- vector<double> obsBuffer;
- vector<double> modelBuffer;
- vector<double> uncertObsBuffer;
+ vector<double> obsLineBuffer;
+ vector<double> obsWindowBuffer;//buffer for observation to calculate average corresponding to model pixel
+ vector<double> model1LineBuffer;
+ vector<double> model2LineBuffer;
+ vector<double> model1buffer;//buffer for model 1 to calculate time regression based on window
+ vector<double> model2buffer;//buffer for model 2 to calculate time regression based on window
+ vector<double> uncertObsLineBuffer;
vector<double> estReadBuffer;
+ vector<double> estWindowBuffer;//buffer for estimate to calculate average corresponding to model pixel
vector<double> uncertReadBuffer;
vector<double> estWriteBuffer(ncol);
vector<double> uncertWriteBuffer(ncol);
for(int irow=0;irow<imgWriterEst.nrOfRow();++irow){
assert(irow<imgReaderEst.nrOfRow());
- imgReaderEst.readData(estReadBuffer,GDT_Float64,irow,0);
+ //not needed here, because we read entire window for each pixel...
+ // imgReaderEst.readData(estReadBuffer,GDT_Float64,irow,0);
imgReaderEst.readData(uncertReadBuffer,GDT_Float64,irow,1);
//read model2 in case current estimate is nodata
imgReaderEst.image2geo(0,irow,x,y);
imgReaderModel2.geo2image(x,y,modCol,modRow);
assert(modRow>=0&&modRow<imgReaderModel2.nrOfRow());
- imgReaderModel2.readData(modelBuffer,GDT_Float64,modRow,0);
+ imgReaderModel2.readData(model2LineBuffer,GDT_Float64,modRow,0);
+
+ imgReaderModel1.geo2image(x,y,modCol,modRow);
+ assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
+ imgReaderModel1.readData(model1LineBuffer,GDT_Float64,modRow,0);
+
if(update){
- imgReaderObs.readData(obsBuffer,GDT_Float64,irow,0);
+ imgReaderObs.readData(obsLineBuffer,GDT_Float64,irow,0);
if(imgReaderObs.nrOfBand()>1)
- imgReaderObs.readData(uncertObsBuffer,GDT_Float64,irow,1);
+ imgReaderObs.readData(uncertObsLineBuffer,GDT_Float64,irow,1);
}
for(int icol=0;icol<imgWriterEst.nrOfCol();++icol){
- double estValue=estReadBuffer[icol];
+ int minCol=(icol>down_opt[0]/2) ? icol-down_opt[0]/2 : 0;
+ int maxCol=(icol+down_opt[0]/2<imgReaderEst.nrOfCol()) ? icol+down_opt[0]/2 : imgReaderEst.nrOfCol()-1;
+ int minRow=(irow>down_opt[0]/2) ? irow-down_opt[0]/2 : 0;
+ int maxRow=(irow+down_opt[0]/2<imgReaderEst.nrOfRow()) ? irow+down_opt[0]/2 : imgReaderEst.nrOfRow()-1;
+ imgReaderEst.readDataBlock(estWindowBuffer,GDT_Float64,minCol,maxCol,minRow,maxRow,0);
+ double estMeanValue=stat.mean(estWindowBuffer);
+ // double estValue=estReadBuffer[icol];
+ double estValue=estWindowBuffer[estWindowBuffer.size()/2];
//time update
- if(imgReaderEst.isNoData(estValue)){
- //pk: in case we have not found any valid data yet, better here to take the current model value
+ imgReaderEst.image2geo(0,irow,x,y);
+ imgReaderModel1.geo2image(x,y,modCol,modRow);
+ double difference=(estMeanValue-c0obs)/c1obs-model1LineBuffer[modCol];
+ difference*=difference;
+ bool estNodata=(sqrt(difference)>uncertModel_opt[0]*stdDev);
+ estNodata=estNodata||imgReaderEst.isNoData(estValue);
+ if(estNodata){
imgReaderEst.image2geo(icol,irow,x,y);
imgReaderModel2.geo2image(x,y,modCol,modRow);
assert(modRow>=0&&modRow<imgReaderModel2.nrOfRow());
- if(imgReaderModel2.isNoData(modelBuffer[modCol])){//if both estimate and model are no-data, set obs to nodata
+ //pk: in case we have not found any valid data yet, better here to take the current model value
+ if(imgReaderModel2.isNoData(model2LineBuffer[modCol])){//if both estimate and model are no-data, set obs to nodata
estWriteBuffer[icol]=obsnodata_opt[0];
uncertWriteBuffer[icol]=uncertNodata_opt[0];
}
else{
- estWriteBuffer[icol]=modelBuffer[modCol];
+ estWriteBuffer[icol]=model2LineBuffer[modCol];
uncertWriteBuffer[icol]=uncertModel_opt[0]*stdDev;
}
}
else{
+ if(window_opt[0]>0){
+ try{
+ imgReaderEst.image2geo(icol,irow,x,y);
+ imgReaderModel1.geo2image(x,y,modCol,modRow);
+ assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
+ minCol=(modCol>window_opt[0]/2) ? modCol-window_opt[0]/2 : 0;
+ maxCol=(modCol+window_opt[0]/2<imgReaderModel1.nrOfCol()) ? modCol+window_opt[0]/2 : imgReaderModel1.nrOfCol()-1;
+ minRow=(modRow>window_opt[0]/2) ? modRow-window_opt[0]/2 : 0;
+ maxRow=(modRow+window_opt[0]/2<imgReaderModel1.nrOfRow()) ? modRow+window_opt[0]/2 : imgReaderModel1.nrOfRow()-1;
+ imgReaderModel1.readDataBlock(model1buffer,GDT_Float64,minCol,maxCol,minRow,maxRow,0);
+ imgReaderEst.image2geo(icol,irow,x,y);
+
+ imgReaderModel2.geo2image(x,y,modCol,modRow);
+ assert(modRow>=0&&modRow<imgReaderModel2.nrOfRow());
+ minCol=(modCol>window_opt[0]/2) ? modCol-window_opt[0]/2 : 0;
+ maxCol=(modCol+window_opt[0]/2<imgReaderModel2.nrOfCol()) ? modCol+window_opt[0]/2 : imgReaderModel2.nrOfCol()-1;
+ minRow=(modRow>window_opt[0]/2) ? modRow-window_opt[0]/2 : 0;
+ maxRow=(modRow+window_opt[0]/2<imgReaderModel2.nrOfRow()) ? modRow+window_opt[0]/2 : imgReaderModel2.nrOfRow()-1;
+ imgReaderModel2.readDataBlock(model2buffer,GDT_Float64,minCol,maxCol,minRow,maxRow,0);
+ }
+ catch(string errorString){
+ cerr << "Error reading data block for " << minCol << "-" << maxCol << ", " << minRow << "-" << maxRow << endl;
+ }
+ //todo: erase non-similar data to observation...
+ vector<double>::iterator it1=model1buffer.begin();
+ vector<double>::iterator it2=model2buffer.begin();
+ while(it1!=model1buffer.end()&&it2!=model2buffer.end()){
+ //erase nodata
+ double difference=(estMeanValue-c0obs)/c1obs-*it1;
+ difference*=difference;
+ bool modNodata=(sqrt(difference)>uncertModel_opt[0]*stdDev);
+ modNodata=modNodata||imgReaderModel1.isNoData(*it1);
+ modNodata=modNodata||imgReaderModel2.isNoData(*it2);
+
+ if(modNodata){
+ model1buffer.erase(it1);
+ model2buffer.erase(it2);
+ }
+ else{
+ ++it1;
+ ++it2;
+ }
+ }
+ if(model1buffer.size()>minreg_opt[0]&&model2buffer.size()>minreg_opt[0])
+ errMod=stat.linear_regression_err(model1buffer,model2buffer,c0mod,c1mod);
+ else{//use global regression...
+ c0mod=c0modGlobal;
+ c1mod=c1modGlobal;
+ }
+ }
double certNorm=(errMod*errMod+errObs*errObs);
double certMod=errObs*errObs/certNorm;
double certObs=errMod*errMod/certNorm;
double regTime=(c0mod+c1mod*estValue)*certMod;
+ //todo: check? estValue is already in the correct scaling from last iteration, see mail to Fer
double regSensor=(c0obs+c1obs*estValue)*certObs;
estWriteBuffer[icol]=regTime+regSensor;
double totalUncertainty=0;
@@ -504,15 +587,15 @@ int main(int argc,char **argv) {
uncertWriteBuffer[icol]=totalUncertainty+uncertReadBuffer[icol];
}
//observation update
- if(update&&!imgReaderObs.isNoData(obsBuffer[icol])){
+ if(update&&!imgReaderObs.isNoData(obsLineBuffer[icol])){
double kalmanGain=1;
double uncertObs=uncertObs_opt[0];
- if(uncertObsBuffer.size()>icol)
- uncertObs=uncertObsBuffer[icol];
+ if(uncertObsLineBuffer.size()>icol)
+ uncertObs=uncertObsLineBuffer[icol];
if((uncertWriteBuffer[icol]+uncertObs)>eps_opt[0])
kalmanGain=uncertWriteBuffer[icol]/(uncertWriteBuffer[icol]+uncertObs);
assert(kalmanGain<=1);
- estWriteBuffer[icol]+=kalmanGain*(obsBuffer[icol]-estWriteBuffer[icol]);
+ estWriteBuffer[icol]+=kalmanGain*(obsLineBuffer[icol]-estWriteBuffer[icol]);
uncertWriteBuffer[icol]*=(1-kalmanGain);
}
}
@@ -589,7 +672,6 @@ int main(int argc,char **argv) {
//write last model as output
if(verbose_opt[0])
cout << "write last model as output" << endl;
- // for(int irow=0;irow<imgWriterEst.nrOfRow();++irow){
for(int irow=0;irow<nrow;++irow){
vector<double> estReadBuffer;
vector<double> estWriteBuffer(ncol);
@@ -635,10 +717,10 @@ int main(int argc,char **argv) {
imgReaderModel1.readData(estReadBuffer,GDT_Float64,modRow);
vector<double> estWriteBuffer(ncol);
vector<double> uncertWriteBuffer(ncol);
- vector<double> uncertObsBuffer;
+ vector<double> uncertObsLineBuffer;
imgReaderObs.readData(estWriteBuffer,GDT_Float64,irow,0);
if(imgReaderObs.nrOfBand()>1)
- imgReaderObs.readData(uncertObsBuffer,GDT_Float64,irow,1);
+ imgReaderObs.readData(uncertObsLineBuffer,GDT_Float64,irow,1);
for(int icol=0;icol<imgWriterEst.nrOfCol();++icol){
if(imgReaderObs.isNoData(estWriteBuffer[icol])){
imgWriterEst.image2geo(icol,irow,x,y);
@@ -652,8 +734,8 @@ int main(int argc,char **argv) {
}
else{
double uncertObs=uncertObs_opt[0];
- if(uncertObsBuffer.size()>icol)
- uncertObs=uncertObsBuffer[icol];
+ if(uncertObsLineBuffer.size()>icol)
+ uncertObs=uncertObsLineBuffer[icol];
if(uncertModel_opt[0]*stdDev+uncertObs>eps_opt[0]){
imgWriterEst.image2geo(icol,irow,x,y);
imgReaderModel1.geo2image(x,y,modCol,modRow);
@@ -722,13 +804,12 @@ int main(int argc,char **argv) {
//to keep it general, we must redo it (overlap might have changed)
pfnProgress(progress,pszMessage,pProgressArg);
- double c0mod=0;
- double c1mod=0;
if(verbose_opt[0])
cout << "Calculating regression for " << imgReaderModel1.getFileName() << " " << imgReaderModel2.getFileName() << endl;
- double errMod=imgreg.getRMSE(imgReaderModel1,imgReaderModel2,c0mod,c1mod);
- // double errMod=imgreg.getRMSE(imgReaderModel1,imgReaderModel2,c0mod,c1mod,verbose_opt[0]);
+
+ double errMod=imgreg.getRMSE(imgReaderModel1,imgReaderModel2,c0modGlobal,c1modGlobal);
+ // double errMod=imgreg.getRMSE(imgReaderModel1,imgReaderModel2,c0modGlobal,c1modGlobal,verbose_opt[0]);
bool update=false;
if(obsindex<relobsindex.size()){
@@ -747,6 +828,8 @@ int main(int argc,char **argv) {
if(verbose_opt[0])
cout << "Calculating regression for " << imgReaderModel2.getFileName() << " " << imgReaderObs.getFileName() << endl;
errObs=imgreg.getRMSE(imgReaderModel2,imgReaderObs,c0obs,c1obs,verbose_opt[0]);
+ if(verbose_opt[0])
+ cout << "c0obs, c1obs: " << c0obs << ", " << c1obs << endl;
}
//prediction (also to fill cloudy pixels in update mode)
string input;
@@ -762,57 +845,132 @@ int main(int argc,char **argv) {
imgReaderEst.setOffset(obsoffset_opt[0]);
imgReaderEst.setScale(obsscale_opt[0]);
- vector<double> obsBuffer;
- vector<double> modelBuffer;
- vector<double> uncertObsBuffer;
+ vector<double> obsLineBuffer;
+ vector<double> obsWindowBuffer;//buffer for observation to calculate average corresponding to model pixel
+ vector<double> model1LineBuffer;
+ vector<double> model2LineBuffer;
+ vector<double> model1buffer;//buffer for model 1 to calculate time regression based on window
+ vector<double> model2buffer;//buffer for model 2 to calculate time regression based on window
+ vector<double> uncertObsLineBuffer;
vector<double> estReadBuffer;
+ vector<double> estWindowBuffer;//buffer for estimate to calculate average corresponding to model pixel
vector<double> uncertReadBuffer;
vector<double> estWriteBuffer(ncol);
vector<double> uncertWriteBuffer(ncol);
-
+
for(int irow=0;irow<imgWriterEst.nrOfRow();++irow){
assert(irow<imgReaderEst.nrOfRow());
- imgReaderEst.readData(estReadBuffer,GDT_Float64,irow,0);
+ //not needed here, because we read entire window for each pixel...
+ // imgReaderEst.readData(estReadBuffer,GDT_Float64,irow,0);
imgReaderEst.readData(uncertReadBuffer,GDT_Float64,irow,1);
//read model2 in case current estimate is nodata
imgReaderEst.image2geo(0,irow,x,y);
imgReaderModel2.geo2image(x,y,modCol,modRow);
assert(modRow>=0&&modRow<imgReaderModel2.nrOfRow());
- imgReaderModel2.readData(modelBuffer,GDT_Float64,modRow,0);
+ imgReaderModel2.readData(model2LineBuffer,GDT_Float64,modRow,0);
+
+ imgReaderModel1.geo2image(x,y,modCol,modRow);
+ assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
+ imgReaderModel1.readData(model1LineBuffer,GDT_Float64,modRow,0);
+
if(update){
- imgReaderObs.readData(obsBuffer,GDT_Float64,irow,0);
+ imgReaderObs.readData(obsLineBuffer,GDT_Float64,irow,0);
if(imgReaderObs.nrOfBand()>1)
- imgReaderObs.readData(uncertObsBuffer,GDT_Float64,irow,1);
+ imgReaderObs.readData(uncertObsLineBuffer,GDT_Float64,irow,1);
}
for(int icol=0;icol<imgWriterEst.nrOfCol();++icol){
- double estValue=estReadBuffer[icol];
+ int minCol=(icol>down_opt[0]/2) ? icol-down_opt[0]/2 : 0;
+ int maxCol=(icol+down_opt[0]/2<imgReaderEst.nrOfCol()) ? icol+down_opt[0]/2 : imgReaderEst.nrOfCol()-1;
+ int minRow=(irow>down_opt[0]/2) ? irow-down_opt[0]/2 : 0;
+ int maxRow=(irow+down_opt[0]/2<imgReaderEst.nrOfRow()) ? irow+down_opt[0]/2 : imgReaderEst.nrOfRow()-1;
+ imgReaderEst.readDataBlock(estWindowBuffer,GDT_Float64,minCol,maxCol,minRow,maxRow,0);
+ double estMeanValue=stat.mean(estWindowBuffer);
+ // double estValue=estReadBuffer[icol];
+ double estValue=estWindowBuffer[estWindowBuffer.size()/2];
//time update
- if(imgReaderEst.isNoData(estValue)){
- //pk: in case we have not found any valid data yet, better here to take the current model value
+ imgReaderEst.image2geo(0,irow,x,y);
+ imgReaderModel1.geo2image(x,y,modCol,modRow);
+ double difference=(estMeanValue-c0obs)/c1obs-model1LineBuffer[modCol];
+ difference*=difference;
+ bool estNodata=(sqrt(difference)>uncertModel_opt[0]*stdDev);
+ estNodata=estNodata||imgReaderEst.isNoData(estValue);
+ if(estNodata){
imgReaderEst.image2geo(icol,irow,x,y);
imgReaderModel2.geo2image(x,y,modCol,modRow);
assert(modRow>=0&&modRow<imgReaderModel2.nrOfRow());
- if(imgReaderModel2.isNoData(modelBuffer[modCol])){//if both estimate and model are no-data, set obs to nodata
+ //pk: in case we have not found any valid data yet, better here to take the current model value
+ if(imgReaderModel2.isNoData(model2LineBuffer[modCol])){//if both estimate and model are no-data, set obs to nodata
estWriteBuffer[icol]=obsnodata_opt[0];
uncertWriteBuffer[icol]=uncertNodata_opt[0];
}
else{
- estWriteBuffer[icol]=modelBuffer[modCol];
+ estWriteBuffer[icol]=model2LineBuffer[modCol];
uncertWriteBuffer[icol]=uncertModel_opt[0]*stdDev;
}
}
else{
+ if(window_opt[0]>0){
+ try{
+ imgReaderEst.image2geo(icol,irow,x,y);
+ imgReaderModel1.geo2image(x,y,modCol,modRow);
+ assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
+ minCol=(modCol>window_opt[0]/2) ? modCol-window_opt[0]/2 : 0;
+ maxCol=(modCol+window_opt[0]/2<imgReaderModel1.nrOfCol()) ? modCol+window_opt[0]/2 : imgReaderModel1.nrOfCol()-1;
+ minRow=(modRow>window_opt[0]/2) ? modRow-window_opt[0]/2 : 0;
+ maxRow=(modRow+window_opt[0]/2<imgReaderModel1.nrOfRow()) ? modRow+window_opt[0]/2 : imgReaderModel1.nrOfRow()-1;
+ imgReaderModel1.readDataBlock(model1buffer,GDT_Float64,minCol,maxCol,minRow,maxRow,0);
+ imgReaderEst.image2geo(icol,irow,x,y);
+
+ imgReaderModel2.geo2image(x,y,modCol,modRow);
+ assert(modRow>=0&&modRow<imgReaderModel2.nrOfRow());
+ minCol=(modCol>window_opt[0]/2) ? modCol-window_opt[0]/2 : 0;
+ maxCol=(modCol+window_opt[0]/2<imgReaderModel2.nrOfCol()) ? modCol+window_opt[0]/2 : imgReaderModel2.nrOfCol()-1;
+ minRow=(modRow>window_opt[0]/2) ? modRow-window_opt[0]/2 : 0;
+ maxRow=(modRow+window_opt[0]/2<imgReaderModel2.nrOfRow()) ? modRow+window_opt[0]/2 : imgReaderModel2.nrOfRow()-1;
+ imgReaderModel2.readDataBlock(model2buffer,GDT_Float64,minCol,maxCol,minRow,maxRow,0);
+ }
+ catch(string errorString){
+ cerr << "Error reading data block for " << minCol << "-" << maxCol << ", " << minRow << "-" << maxRow << endl;
+ }
+ //todo: erase non-similar data to observation...
+ vector<double>::iterator it1=model1buffer.begin();
+ vector<double>::iterator it2=model2buffer.begin();
+ while(it1!=model1buffer.end()&&it2!=model2buffer.end()){
+ //erase nodata
+ double difference=(estMeanValue-c0obs)/c1obs-*it1;
+ difference*=difference;
+ bool modNodata=(sqrt(difference)>uncertModel_opt[0]*stdDev);
+ modNodata=modNodata||imgReaderModel1.isNoData(*it1);
+ modNodata=modNodata||imgReaderModel2.isNoData(*it2);
+
+ if(modNodata){
+ model1buffer.erase(it1);
+ model2buffer.erase(it2);
+ }
+ else{
+ ++it1;
+ ++it2;
+ }
+ }
+ if(model1buffer.size()>minreg_opt[0]&&model2buffer.size()>minreg_opt[5])
+ errMod=stat.linear_regression_err(model1buffer,model2buffer,c0mod,c1mod);
+ else{//use global regression...
+ c0mod=c0modGlobal;
+ c1mod=c1modGlobal;
+ }
+ }
double certNorm=(errMod*errMod+errObs*errObs);
double certMod=errObs*errObs/certNorm;
double certObs=errMod*errMod/certNorm;
double regTime=(c0mod+c1mod*estValue)*certMod;
+ //todo: check? estValue is already in the correct scaling from last iteration, see mail to Fer
double regSensor=(c0obs+c1obs*estValue)*certObs;
estWriteBuffer[icol]=regTime+regSensor;
double totalUncertainty=0;
if(errMod<eps_opt[0])
totalUncertainty=errObs;
else if(errObs<eps_opt[0])
- totalUncertainty=errObs;
+ totalUncertainty=errMod;
else{
totalUncertainty=1.0/errMod/errMod+1/errObs/errObs;
totalUncertainty=sqrt(1.0/totalUncertainty);
@@ -820,15 +978,15 @@ int main(int argc,char **argv) {
uncertWriteBuffer[icol]=totalUncertainty+uncertReadBuffer[icol];
}
//observation update
- if(update&&!imgReaderObs.isNoData(obsBuffer[icol])){
+ if(update&&!imgReaderObs.isNoData(obsLineBuffer[icol])){
double kalmanGain=1;
double uncertObs=uncertObs_opt[0];
- if(uncertObsBuffer.size()>icol)
- uncertObs=uncertObsBuffer[icol];
+ if(uncertObsLineBuffer.size()>icol)
+ uncertObs=uncertObsLineBuffer[icol];
if((uncertWriteBuffer[icol]+uncertObs)>eps_opt[0])
kalmanGain=uncertWriteBuffer[icol]/(uncertWriteBuffer[icol]+uncertObs);
assert(kalmanGain<=1);
- estWriteBuffer[icol]+=kalmanGain*(obsBuffer[icol]-estWriteBuffer[icol]);
+ estWriteBuffer[icol]+=kalmanGain*(obsLineBuffer[icol]-estWriteBuffer[icol]);
uncertWriteBuffer[icol]*=(1-kalmanGain);
}
}
@@ -907,7 +1065,7 @@ int main(int argc,char **argv) {
vector<double> estForwardBuffer;
vector<double> estBackwardBuffer;
- vector<double> uncertObsBuffer;
+ vector<double> uncertObsLineBuffer;
vector<double> uncertForwardBuffer;
vector<double> uncertBackwardBuffer;
vector<double> uncertReadBuffer;
@@ -943,7 +1101,7 @@ int main(int argc,char **argv) {
if(update){
imgReaderObs.readData(estWriteBuffer,GDT_Float64,irow,0);
if(imgReaderObs.nrOfBand()>1)
- imgReaderObs.readData(uncertObsBuffer,GDT_Float64,irow,1);
+ imgReaderObs.readData(uncertObsLineBuffer,GDT_Float64,irow,1);
}
for(int icol=0;icol<imgWriterEst.nrOfCol();++icol){
@@ -956,8 +1114,8 @@ int main(int argc,char **argv) {
// if(update){//check for nodata in observation
// if(imgReaderObs.isNoData(estWriteBuffer[icol]))
// uncertObs=uncertNodata_opt[0];
- // else if(uncertObsBuffer.size()>icol)
- // uncertObs=uncertObsBuffer[icol];
+ // else if(uncertObsLineBuffer.size()>icol)
+ // uncertObs=uncertObsLineBuffer[icol];
// }
double noemer=(C+D);
diff --git a/src/imageclasses/ImgReaderGdal.h b/src/imageclasses/ImgReaderGdal.h
index 5d4fceb..d287631 100644
--- a/src/imageclasses/ImgReaderGdal.h
+++ b/src/imageclasses/ImgReaderGdal.h
@@ -237,14 +237,25 @@ template<typename T> void ImgReaderGdal::readDataBlock(std::vector<T>& buffer, c
GDALRasterBand *poBand;
assert(band<nrOfBand()+1);
poBand = m_gds->GetRasterBand(band+1);//GDAL uses 1 based index
- assert(minCol<nrOfCol());
- assert(minCol>=0);
- assert(maxCol<nrOfCol());
- assert(minCol<=maxCol);
- assert(minRow<nrOfRow());
- assert(minRow>=0);
- assert(maxRow<nrOfRow());
- assert(minRow<=maxRow);
+ if(minCol>=nrOfCol() ||
+ (minCol<0) ||
+ (maxCol>=nrOfCol()) ||
+ (minCol>maxCol) ||
+ (minRow>=nrOfRow()) ||
+ (minRow<0) ||
+ (maxRow>=nrOfRow()) ||
+ (minRow>maxRow)){
+ std::string errorString="block not within image boundaries";
+ throw(errorString);
+ }
+ /* assert(minCol<nrOfCol()); */
+ /* assert(minCol>=0); */
+ /* assert(maxCol<nrOfCol()); */
+ /* assert(minCol<=maxCol); */
+ /* assert(minRow<nrOfRow()); */
+ /* assert(minRow>=0); */
+ /* assert(maxRow<nrOfRow()); */
+ /* assert(minRow<=maxRow); */
if(buffer.size()!=(maxRow-minRow+1)*(maxCol-minCol+1))
buffer.resize((maxRow-minRow+1)*(maxCol-minCol+1));
poBand->RasterIO(GF_Read,minCol,minRow,maxCol-minCol+1,maxRow-minRow+1,&(buffer[0]),(maxCol-minCol+1),(maxRow-minRow+1),dataType,0,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