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 55c8d93fbbaf8f5d579c32075ea392523e3f41d8 Author: Pieter Kempeneers <kempe...@gmail.com> Date: Wed Oct 15 18:33:23 2014 +0200 sensor regression based on model data input --- src/apps/pkkalman.cc | 242 ++++++++++++++++++++++++++++++++------------------- 1 file changed, 154 insertions(+), 88 deletions(-) diff --git a/src/apps/pkkalman.cc b/src/apps/pkkalman.cc index 04daa31..bbc7475 100644 --- a/src/apps/pkkalman.cc +++ b/src/apps/pkkalman.cc @@ -56,6 +56,8 @@ int main(int argc,char **argv) { 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<bool> regObs_opt("regobs", "regobs", "Perform regression between modeled and observed value",false); + Optionpk<bool> checkDiff_opt("diff", "diff", "Flag observation as invalid if difference with model is above uncertainty",false); 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."); @@ -87,6 +89,8 @@ int main(int argc,char **argv) { down_opt.retrieveOption(argc,argv); threshold_opt.retrieveOption(argc,argv); minreg_opt.retrieveOption(argc,argv); + // regObs_opt.retrieveOption(argc,argv); + checkDiff_opt.retrieveOption(argc,argv); window_opt.retrieveOption(argc,argv); oformat_opt.retrieveOption(argc,argv); option_opt.retrieveOption(argc,argv); @@ -462,7 +466,7 @@ int main(int argc,char **argv) { 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> estWindowBuffer;//buffer for estimate to calculate average corresponding to model pixel vector<double> uncertReadBuffer; vector<double> estWriteBuffer(ncol); vector<double> uncertWriteBuffer(ncol); @@ -470,7 +474,7 @@ int main(int argc,char **argv) { for(int irow=0;irow<imgWriterEst.nrOfRow();++irow){ assert(irow<imgReaderEst.nrOfRow()); //not needed here, because we read entire window for each pixel... - // imgReaderEst.readData(estReadBuffer,GDT_Float64,irow,0); + 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); @@ -492,35 +496,55 @@ int main(int argc,char **argv) { 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]; + if(update) + imgReaderObs.readDataBlock(obsWindowBuffer,GDT_Float64,minCol,maxCol,minRow,maxRow,0); + // imgReaderEst.readDataBlock(estWindowBuffer,GDT_Float64,minCol,maxCol,minRow,maxRow,0); + double estValue=estReadBuffer[icol]; + // double estValue=estWindowBuffer[estWindowBuffer.size()/2]; + double estMeanValue=0;//stat.mean(estWindowBuffer); + double nvalid=0; //time update - 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); + imgReaderEst.image2geo(icol,irow,x,y); + imgReaderModel2.geo2image(x,y,modCol,modRow); + assert(modRow>=0&&modRow<imgReaderModel2.nrOfRow()); + double modValue=model2LineBuffer[modCol]; + bool estNodata=imgReaderEst.isNoData(estValue); + //todo: debug checkDiff_opt + // if(checkDiff_opt[0]){ + // vector<double>::iterator itwin=estWindowBuffer.begin(); + // while(itwin!=estWindowBuffer.end()){ + // if(!imgReaderEst.isNoData(*itwin)){ + // estMeanValue+=*itwin; + // ++nvalid; + // } + // ++itwin; + // } + // estMeanValue/=nvalid; + // double difference=(estMeanValue-c0obs)/c1obs-model1LineBuffer[modCol]; + // difference*=difference; + // estNodata=estNodata||(sqrt(difference)>uncertModel_opt[0]*stdDev); + // } if(estNodata){ - imgReaderEst.image2geo(icol,irow,x,y); - imgReaderModel2.geo2image(x,y,modCol,modRow); - assert(modRow>=0&&modRow<imgReaderModel2.nrOfRow()); - //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 + //we have not found any valid data yet, better here to take the current model value if valid + if(imgReaderModel2.isNoData(modValue)){//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]=model2LineBuffer[modCol]; + estWriteBuffer[icol]=modValue; uncertWriteBuffer[icol]=uncertModel_opt[0]*stdDev; } } else{ if(window_opt[0]>0){ try{ - imgReaderEst.image2geo(icol,irow,x,y); + // imgReaderModel2.geo2image(x,y,modCol,modRow);//did that already + 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); + imgReaderModel1.geo2image(x,y,modCol,modRow); assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow()); minCol=(modCol>window_opt[0]/2) ? modCol-window_opt[0]/2 : 0; @@ -528,30 +552,25 @@ int main(int argc,char **argv) { 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); + // imgReaderEst.image2geo(icol,irow,x,y); } catch(string errorString){ cerr << "Error reading data block for " << minCol << "-" << maxCol << ", " << minRow << "-" << maxRow << endl; } - //todo: erase non-similar data to observation... + //erase no-data from buffer 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); + bool modNodata=false; modNodata=modNodata||imgReaderModel1.isNoData(*it1); modNodata=modNodata||imgReaderModel2.isNoData(*it2); - + //todo: debug checkDiff_opt + // if(checkDiff_opt[0]){ + // double difference=(estMeanValue-c0obs)/c1obs-*it1; + // difference*=difference; + // modNodata=modNodata||(sqrt(difference)>uncertModel_opt[0]*stdDev); + // } if(modNodata){ model1buffer.erase(it1); model2buffer.erase(it2); @@ -568,12 +587,17 @@ int main(int argc,char **argv) { c1mod=c1modGlobal; } } + else{ + 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; + + // double regSensor=(c0obs+c1obs*estValue)*certObs; + double regSensor=(c0obs+c1obs*modValue)*certObs; estWriteBuffer[icol]=regTime+regSensor; double totalUncertainty=0; if(errMod<eps_opt[0]) @@ -588,15 +612,26 @@ int main(int argc,char **argv) { } //observation update if(update&&!imgReaderObs.isNoData(obsLineBuffer[icol])){ - double kalmanGain=1; - double uncertObs=uncertObs_opt[0]; - 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*(obsLineBuffer[icol]-estWriteBuffer[icol]); - uncertWriteBuffer[icol]*=(1-kalmanGain); + bool doUpdate=true; + if(checkDiff_opt[0]){ + statfactory::StatFactory statobs; + statobs.setNoDataValues(obsnodata_opt); + double obsMeanValue=statobs.mean(obsWindowBuffer); + double difference=(obsMeanValue-c0obs)/c1obs-modValue; + difference*=difference; + doUpdate=(sqrt(difference)<uncertModel_opt[0]*stdDev); + } + if(doUpdate){ + double kalmanGain=1; + double uncertObs=uncertObs_opt[0]; + 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*(obsLineBuffer[icol]-estWriteBuffer[icol]); + uncertWriteBuffer[icol]*=(1-kalmanGain); + } } } imgWriterEst.writeData(estWriteBuffer,GDT_Float64,irow,0); @@ -853,7 +888,7 @@ int main(int argc,char **argv) { 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> estWindowBuffer;//buffer for estimate to calculate average corresponding to model pixel vector<double> uncertReadBuffer; vector<double> estWriteBuffer(ncol); vector<double> uncertWriteBuffer(ncol); @@ -861,7 +896,7 @@ int main(int argc,char **argv) { for(int irow=0;irow<imgWriterEst.nrOfRow();++irow){ assert(irow<imgReaderEst.nrOfRow()); //not needed here, because we read entire window for each pixel... - // imgReaderEst.readData(estReadBuffer,GDT_Float64,irow,0); + 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); @@ -883,35 +918,55 @@ int main(int argc,char **argv) { 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]; + if(update) + imgReaderObs.readDataBlock(obsWindowBuffer,GDT_Float64,minCol,maxCol,minRow,maxRow,0); + // imgReaderEst.readDataBlock(estWindowBuffer,GDT_Float64,minCol,maxCol,minRow,maxRow,0); + double estValue=estReadBuffer[icol]; + // double estValue=estWindowBuffer[estWindowBuffer.size()/2]; + double estMeanValue=0;//stat.mean(estWindowBuffer); + double nvalid=0; //time update - 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); + imgReaderEst.image2geo(icol,irow,x,y); + imgReaderModel2.geo2image(x,y,modCol,modRow); + assert(modRow>=0&&modRow<imgReaderModel2.nrOfRow()); + double modValue=model2LineBuffer[modCol]; + bool estNodata=imgReaderEst.isNoData(estValue); + //todo: debug checkDiff_opt + // if(checkDiff_opt[0]){ + // vector<double>::iterator itwin=estWindowBuffer.begin(); + // while(itwin!=estWindowBuffer.end()){ + // if(!imgReaderEst.isNoData(*itwin)){ + // estMeanValue+=*itwin; + // ++nvalid; + // } + // ++itwin; + // } + // estMeanValue/=nvalid; + // double difference=(estMeanValue-c0obs)/c1obs-model1LineBuffer[modCol]; + // difference*=difference; + // estNodata=estNodata||(sqrt(difference)>uncertModel_opt[0]*stdDev); + // } if(estNodata){ - imgReaderEst.image2geo(icol,irow,x,y); - imgReaderModel2.geo2image(x,y,modCol,modRow); - assert(modRow>=0&&modRow<imgReaderModel2.nrOfRow()); - //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 + //we have not found any valid data yet, better here to take the current model value if valid + if(imgReaderModel2.isNoData(modValue)){//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]=model2LineBuffer[modCol]; + estWriteBuffer[icol]=modValue; uncertWriteBuffer[icol]=uncertModel_opt[0]*stdDev; } } else{ if(window_opt[0]>0){ try{ - imgReaderEst.image2geo(icol,irow,x,y); + // imgReaderModel2.geo2image(x,y,modCol,modRow);//did that already + 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); + imgReaderModel1.geo2image(x,y,modCol,modRow); assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow()); minCol=(modCol>window_opt[0]/2) ? modCol-window_opt[0]/2 : 0; @@ -919,30 +974,25 @@ int main(int argc,char **argv) { 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); + // imgReaderEst.image2geo(icol,irow,x,y); } catch(string errorString){ cerr << "Error reading data block for " << minCol << "-" << maxCol << ", " << minRow << "-" << maxRow << endl; } - //todo: erase non-similar data to observation... + //erase no-data from buffer 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); + bool modNodata=false; modNodata=modNodata||imgReaderModel1.isNoData(*it1); modNodata=modNodata||imgReaderModel2.isNoData(*it2); - + //todo: debug checkDiff_opt + // if(checkDiff_opt[0]){ + // double difference=(estMeanValue-c0obs)/c1obs-*it1; + // difference*=difference; + // modNodata=modNodata||(sqrt(difference)>uncertModel_opt[0]*stdDev); + // } if(modNodata){ model1buffer.erase(it1); model2buffer.erase(it2); @@ -959,12 +1009,17 @@ int main(int argc,char **argv) { c1mod=c1modGlobal; } } + else{ + 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; + + // double regSensor=(c0obs+c1obs*estValue)*certObs; + double regSensor=(c0obs+c1obs*modValue)*certObs; estWriteBuffer[icol]=regTime+regSensor; double totalUncertainty=0; if(errMod<eps_opt[0]) @@ -979,15 +1034,26 @@ int main(int argc,char **argv) { } //observation update if(update&&!imgReaderObs.isNoData(obsLineBuffer[icol])){ - double kalmanGain=1; - double uncertObs=uncertObs_opt[0]; - 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*(obsLineBuffer[icol]-estWriteBuffer[icol]); - uncertWriteBuffer[icol]*=(1-kalmanGain); + bool doUpdate=true; + if(checkDiff_opt[0]){ + statfactory::StatFactory statobs; + statobs.setNoDataValues(obsnodata_opt); + double obsMeanValue=statobs.mean(obsWindowBuffer); + double difference=(obsMeanValue-c0obs)/c1obs-modValue; + difference*=difference; + doUpdate=(sqrt(difference)<uncertModel_opt[0]*stdDev); + } + if(doUpdate){ + double kalmanGain=1; + double uncertObs=uncertObs_opt[0]; + 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*(obsLineBuffer[icol]-estWriteBuffer[icol]); + uncertWriteBuffer[icol]*=(1-kalmanGain); + } } } imgWriterEst.writeData(estWriteBuffer,GDT_Float64,irow,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