[pktools] 181/375: hist in StatFactory
Bas Couwenberg
sebastic at xs4all.nl
Wed Dec 3 21:54:12 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 033dfe059aba250ab09740ef87b0c891a83ab0ad
Author: Pieter Kempeneers <kempenep at gmail.com>
Date: Sun Jan 12 14:34:27 2014 +0100
hist in StatFactory
---
src/algorithms/StatFactory.h | 8 ++--
src/apps/pkfilter.cc | 5 +++
src/apps/pkfilterdem.cc | 33 +++++++++-----
src/apps/pkinfo.cc | 2 +
src/apps/pkstatascii.cc | 95 ++++++++++++++++++++++++++++++++-------
src/imageclasses/ImgReaderGdal.cc | 2 +-
6 files changed, 114 insertions(+), 31 deletions(-)
diff --git a/src/algorithms/StatFactory.h b/src/algorithms/StatFactory.h
index 0b95b76..e59b192 100644
--- a/src/algorithms/StatFactory.h
+++ b/src/algorithms/StatFactory.h
@@ -120,17 +120,17 @@ public:
m_noDataValues=vnodata;
return m_noDataValues.size();
};
- double getRandomValue(const gsl_rng* r, const std::string type, double a=0, double b=0) const{
+ double getRandomValue(const gsl_rng* r, const std::string type, double a=0, double b=1) const{
std::map<std::string, DISTRIBUTION_TYPE> m_distMap;
initDist(m_distMap);
double randValue=0;
switch(m_distMap[type]){
case(uniform):
- randValue = gsl_rng_uniform(r);
+ randValue = a+gsl_rng_uniform(r);
break;
case(gaussian):
default:
- randValue = gsl_ran_gaussian(r, a);
+ randValue = a+gsl_ran_gaussian(r, b);
break;
}
return randValue;
@@ -619,7 +619,7 @@ template<class T> void StatFactory::distribution(const std::vector<T>& input, t
if(*it==maximum)
theBin=nbin-1;
else if(*it>minimum && *it<maximum)
- theBin=static_cast<int>(static_cast<double>((*it)-minimum)/(maximum-minimum)*nbin);
+ theBin=static_cast<int>(static_cast<double>((nbin-1)*(*it)-minimum)/(maximum-minimum));
++output[theBin];
// if(*it==maximum)
// ++output[nbin-1];
diff --git a/src/apps/pkfilter.cc b/src/apps/pkfilter.cc
index 0efcd31..82666fb 100644
--- a/src/apps/pkfilter.cc
+++ b/src/apps/pkfilter.cc
@@ -210,6 +210,11 @@ int main(int argc,char **argv) {
}
else if(input.getColorTable()!=NULL)
output.setColorTable(input.getColorTable());
+
+ if(nodata_opt.size()){
+ for(int iband=0;iband<output.nrOfBand();++iband)
+ output.GDALSetNoDataValue(nodata_opt[0],iband);
+ }
filter2d::Filter2d filter2d;
filter::Filter filter1d;
diff --git a/src/apps/pkfilterdem.cc b/src/apps/pkfilterdem.cc
index 77f857e..034d53e 100644
--- a/src/apps/pkfilterdem.cc
+++ b/src/apps/pkfilterdem.cc
@@ -35,7 +35,7 @@ int main(int argc,char **argv) {
Optionpk<std::string> output_opt("o", "output", "Output image file");
Optionpk<std::string> tmpdir_opt("tmp", "tmp", "Temporary directory","/tmp",2);
Optionpk<bool> disc_opt("circ", "circular", "circular disc kernel for dilation and erosion", false);
- Optionpk<string> postFilter_opt("pf", "pfilter", "post processing filter: etew_min, promorph (progressive morphological filter),open,close,none).");
+ Optionpk<string> postFilter_opt("f", "filter", "post processing filter: etew_min, promorph (progressive morphological filter),open,close).");
Optionpk<double> dim_opt("dim", "dim", "maximum filter kernel size (optionally you can set both initial and maximum filter kernel size", 17);
Optionpk<double> maxSlope_opt("st", "st", "slope threshold used for morphological filtering. Use a low values to remove more height objects in flat terrains", 0.0);
Optionpk<double> hThreshold_opt("ht", "ht", "initial height threshold for progressive morphological filtering. Use low values to remove more height objects. Optionally, a maximum height threshold can be set via a second argument (e.g., -ht 0.2 -ht 2.5 sets an initial threshold at 0.2 m and caps the threshold at 2.5 m).", 0.2);
@@ -44,6 +44,7 @@ int main(int argc,char **argv) {
Optionpk<string> oformat_opt("of", "oformat", "Output image format (see also gdal_translate). Empty string: inherit from input image");
Optionpk<string> colorTable_opt("ct", "ct", "color table (file with 5 columns: id R G B ALFA (0: transparent, 255: solid). Use none to ommit color table");
Optionpk<string> option_opt("co", "co", "Creation option for output file. Multiple options can be specified.");
+ Optionpk<short> nodata_opt("nodata", "nodata", "nodata value(s) for smoothnodata filter");
Optionpk<short> verbose_opt("v", "verbose", "verbose mode if > 0", 0);
bool doProcess;//stop process when program was invoked with help option (-h --help)
@@ -60,6 +61,7 @@ int main(int argc,char **argv) {
otype_opt.retrieveOption(argc,argv);
oformat_opt.retrieveOption(argc,argv);
colorTable_opt.retrieveOption(argc,argv);
+ nodata_opt.retrieveOption(argc,argv);
verbose_opt.retrieveOption(argc,argv);
}
catch(string predefinedString){
@@ -122,6 +124,12 @@ int main(int argc,char **argv) {
if(colorTable_opt.size())
outputWriter.setColorTable(colorTable_opt[0]);
+ //set nodata value
+ if(nodata_opt.size()){
+ for(int iband=0;iband<outputWriter.nrOfBand();++iband)
+ outputWriter.GDALSetNoDataValue(nodata_opt[0],iband);
+ }
+
Vector2d<double> inputData(input.nrOfRow(),input.nrOfCol());
Vector2d<double> outputData(outputWriter.nrOfRow(),outputWriter.nrOfCol());
Vector2d<double> tmpData(outputWriter.nrOfRow(),outputWriter.nrOfCol());
@@ -143,18 +151,24 @@ int main(int argc,char **argv) {
dim_opt.insert(dim_opt.begin(),dim_opt[1]);
dim_opt.erase(dim_opt.begin()+2);
}
+
+ filter2d::Filter2d theFilter;
+ if(nodata_opt.size()){
+ for(int inodata=0;inodata<nodata_opt.size();++inodata)
+ theFilter.pushNoDataValue(nodata_opt[inodata]);
+ }
+
unsigned long int nchange=1;
if(postFilter_opt[0]=="etew_min"){
//Elevation Threshold with Expand Window (ETEW) Filter (p.73 from Airborne LIDAR Data Processing and Analysis Tools ALDPAT 1.0)
//first iteration is performed assuming only minima are selected using options -fir all -comp min
//increase cells and thresholds until no points from the previous iteration are discarded.
int dim=dim_opt[0];
- filter2d::Filter2d morphFilter;
- // morphFilter.setNoValue(0);
+ // theFilter.setNoValue(0);
int iteration=1;
while(nchange>minChange_opt[0]&&dim<=dim_opt[1]){
double hThreshold=maxSlope_opt[0]*dim;
- nchange=morphFilter.morphology(inputData,outputData,"erode",dim,dim,disc_opt[0],hThreshold);
+ nchange=theFilter.morphology(inputData,outputData,"erode",dim,dim,disc_opt[0],hThreshold);
inputData=outputData;
dim+=2;//change from theory: originally double cellCize
std::cout << "iteration " << iteration << ": " << nchange << " pixels changed" << std::endl;
@@ -166,7 +180,6 @@ int main(int argc,char **argv) {
//first iteration is performed assuming only minima are selected using options -fir all -comp min
//increase cells and thresholds until no points from the previous iteration are discarded.
int dim=dim_opt[0];
- filter2d::Filter2d theFilter;
double hThreshold=hThreshold_opt[0];
int iteration=1;
while(nchange>minChange_opt[0]&&dim<=dim_opt[1]){
@@ -193,10 +206,9 @@ int main(int argc,char **argv) {
}
}
else if(postFilter_opt[0]=="open"){
- filter2d::Filter2d morphFilter;
try{
- morphFilter.morphology(inputData,tmpData,"erode",dim_opt[0],dim_opt[0],disc_opt[0],maxSlope_opt[0]);
- morphFilter.morphology(tmpData,outputData,"dilate",dim_opt[0],dim_opt[0],disc_opt[0],maxSlope_opt[0]);
+ theFilter.morphology(inputData,tmpData,"erode",dim_opt[0],dim_opt[0],disc_opt[0],hThreshold_opt[0]);
+ theFilter.morphology(tmpData,outputData,"dilate",dim_opt[0],dim_opt[0],disc_opt[0],hThreshold_opt[0]);
outputData=inputData;
}
catch(std::string errorString){
@@ -205,10 +217,9 @@ int main(int argc,char **argv) {
}
}
else if(postFilter_opt[0]=="close"){
- filter2d::Filter2d morphFilter;
try{
- morphFilter.morphology(inputData,tmpData,"dilate",dim_opt[0],dim_opt[0],disc_opt[0],maxSlope_opt[0]);
- morphFilter.morphology(tmpData,outputData,"erode",dim_opt[0],dim_opt[0],disc_opt[0],maxSlope_opt[0]);
+ theFilter.morphology(inputData,tmpData,"dilate",dim_opt[0],dim_opt[0],disc_opt[0],hThreshold_opt[0]);
+ theFilter.morphology(tmpData,outputData,"erode",dim_opt[0],dim_opt[0],disc_opt[0],hThreshold_opt[0]);
}
catch(std::string errorString){
cout << errorString << endl;
diff --git a/src/apps/pkinfo.cc b/src/apps/pkinfo.cc
index a0ffa5a..e2d9294 100644
--- a/src/apps/pkinfo.cc
+++ b/src/apps/pkinfo.cc
@@ -24,6 +24,8 @@ along with pktools. If not, see <http://www.gnu.org/licenses/>.
#include "imageclasses/ImgReaderGdal.h"
#include "imageclasses/ImgReaderOgr.h"
+using namespace std;
+
int main(int argc, char *argv[])
{
Optionpk<std::string> input_opt("i","input","Input image file");
diff --git a/src/apps/pkstatascii.cc b/src/apps/pkstatascii.cc
index 71ba3c2..4cf297d 100644
--- a/src/apps/pkstatascii.cc
+++ b/src/apps/pkstatascii.cc
@@ -28,7 +28,7 @@ using namespace std;
int main(int argc, char *argv[])
{
- Optionpk<string> input_opt("i","input","name of the input text file","");
+ Optionpk<string> input_opt("i","input","name of the input text file");
Optionpk<char> fs_opt("fs","fs","field separator.",' ');
Optionpk<char> comment_opt("comment","comment","comment character",'#');
Optionpk<bool> output_opt("o","output","output the selected columns",false);
@@ -38,8 +38,8 @@ int main(int argc, char *argv[])
Optionpk<bool> size_opt("size","size","sample size",false);
Optionpk<unsigned int> rand_opt("rnd", "rnd", "generate random numbers", 0);
Optionpk<std::string> randdist_opt("dist", "dist", "distribution for generating random numbers, see http://www.gn/software/gsl/manual/gsl-ref_toc.html#TOC320 (only uniform and Gaussian supported yet)", "gaussian");
- Optionpk<double> randa_opt("rnda", "rnda", "first parameter for random distribution (standard deviation in case of Gaussian)", 1);
- Optionpk<double> randb_opt("rndb", "rndb", "second parameter for random distribution (standard deviation in case of Gaussian)", 0);
+ Optionpk<double> randa_opt("rnda", "rnda", "first parameter for random distribution (mean value in case of Gaussian)", 0);
+ Optionpk<double> randb_opt("rndb", "rndb", "second parameter for random distribution (standard deviation in case of Gaussian)", 1);
Optionpk<bool> mean_opt("mean","mean","calculate median",false);
Optionpk<bool> median_opt("median","median","calculate median",false);
Optionpk<bool> var_opt("var","var","calculate variance",false);
@@ -131,28 +131,38 @@ int main(int argc, char *argv[])
assert(dataVector.size());
double minValue=0;
double maxValue=0;
-
- unsigned int nbin=(nbin_opt.size())? nbin_opt[0]:0;
+ unsigned int nbin=0;
+ if(nbin_opt.size())
+ nbin=nbin_opt[0];
if(histogram_opt[0]){
+ stat.minmax(dataVector[0],dataVector[0].begin(),dataVector[0].end(),minValue,maxValue);
+ if(src_min_opt.size())
+ minValue=src_min_opt[0];
+ if(src_max_opt.size())
+ maxValue=src_max_opt[0];
if(nbin<1){
std::cerr << "Warning: number of bins not defined, calculating bins from min and max value" << std::endl;
- if(maxValue<=minValue)
- stat.minmax(dataVector[0],dataVector[0].begin(),dataVector[0].end(),minValue,maxValue);
nbin=maxValue-minValue+1;
}
}
- double minX=src_min_opt[0];
- double minY=(src_min_opt.size()==2)? src_min_opt[1] : src_min_opt[0];
- double maxX=src_max_opt[0];
- double maxY=(src_max_opt.size()==2)? src_max_opt[1] : src_max_opt[0];
+ double minX=0;
+ double minY=0;
+ double maxX=0;
+ double maxY=0;
if(histogram2d_opt[0]){
assert(col_opt.size()==2);
if(nbin<1){
std::cerr << "Warning: number of bins not defined, calculating bins from min and max value" << std::endl;
- if(maxValue<=minValue){
- stat.minmax(dataVector[0],dataVector[0].begin(),dataVector[0].end(),minX,maxX);
- stat.minmax(dataVector[1],dataVector[1].begin(),dataVector[1].end(),minY,maxY);
- }
+ stat.minmax(dataVector[0],dataVector[0].begin(),dataVector[0].end(),minX,maxX);
+ stat.minmax(dataVector[1],dataVector[1].begin(),dataVector[1].end(),minY,maxY);
+ if(src_min_opt.size())
+ minX=src_min_opt[0];
+ if(src_min_opt.size()>1)
+ minY=src_min_opt[1];
+ if(src_max_opt.size())
+ maxX=src_max_opt[0];
+ if(src_max_opt.size()>1)
+ maxY=src_max_opt[1];
minValue=(minX<minY)? minX:minY;
maxValue=(maxX>maxY)? maxX:maxY;
if(verbose_opt[0])
@@ -207,7 +217,62 @@ int main(int argc, char *argv[])
else
std::cout << "calculating histogram for col " << icol << std::endl;
}
+ //test
+ // cout << "debug0" << endl;
+ // cout << "dataVector.size(): " << dataVector.size() << endl;
+ // cout << "statVector.size(): " << statVector.size() << endl;
+
+ // double theMinValue=0;
+ // double theMaxValue=0;
+
+ // stat.minmax(dataVector[icol],dataVector[icol].begin(),dataVector[icol].end(),theMinValue,theMaxValue);
+ // if(minValue<maxValue&&minValue>theMinValue)
+ // theMinValue=minValue;
+ // if(minValue<maxValue&&maxValue<theMaxValue)
+ // theMaxValue=maxValue;
+
+ // //todo: check...
+ // minValue=theMinValue;
+ // maxValue=theMaxValue;
+
+ // if(maxValue<=minValue){
+ // std::ostringstream s;
+ // s<<"Error: could not calculate distribution (min>=max)";
+ // throw(s.str());
+ // }
+ // assert(nbin);
+ // assert(dataVector[icol].size());
+ // if(statVector[icol].size()!=nbin){
+ // statVector[icol].resize(nbin);
+ // for(int i=0;i<nbin;statVector[icol][i++]=0);
+ // }
+ // typename std::vector<double>::const_iterator it;
+ // for(it=dataVector[icol].begin();it!=dataVector[icol].end();++it){
+ // if(*it<minValue)
+ // continue;
+ // if(*it>maxValue)
+ // continue;
+ // if(stat.isNoData(*it))
+ // continue;
+ // int theBin=0;
+ // if(*it==maxValue)
+ // theBin=nbin-1;
+ // else if(*it>minValue && *it<maxValue)
+ // theBin=static_cast<int>(static_cast<double>((nbin-1)*(*it)-minValue)/(maxValue-minValue));
+ // assert(theBin<statVector[icol].size());
+ // ++statVector[icol][theBin];
+ // // if(*it==maxValue)
+ // // ++statVector[icol][nbin-1];
+ // // else if(*it>=minValue && *it<maxValue)
+ // // ++statVector[icol][static_cast<int>(static_cast<double>((*it)-minValue)/(maxValue-minValue)*nbin)];
+ // }
+
+ // exit(0);
+ //end test
+
stat.distribution(dataVector[icol],dataVector[icol].begin(),dataVector[icol].end(),statVector[icol],nbin,minValue,maxValue,sigma);
+ //test
+ cout << "debug1" << endl;
if(verbose_opt[0])
std::cout << "min and max values: " << minValue << ", " << maxValue << std::endl;
}
diff --git a/src/imageclasses/ImgReaderGdal.cc b/src/imageclasses/ImgReaderGdal.cc
index 0c8acc1..7a7cc76 100644
--- a/src/imageclasses/ImgReaderGdal.cc
+++ b/src/imageclasses/ImgReaderGdal.cc
@@ -520,7 +520,7 @@ unsigned long int ImgReaderGdal::getHistogram(std::vector<unsigned long int>& hi
if(maxValue>minValue){
if(nbin==0)
nbin=maxValue-minValue+1;
- scale=static_cast<double>(nbin)/(maxValue-minValue);
+ scale=static_cast<double>(nbin-1)/(maxValue-minValue);
}
else
nbin=1;
--
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