[pktools] 180/375: moved post filter from pklas2img to pkfilterdem
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 b35ec98d112f1d0924b32e207f2ffe0e610eb262
Author: Pieter Kempeneers <kempenep at gmail.com>
Date: Sat Jan 11 23:23:21 2014 +0100
moved post filter from pklas2img to pkfilterdem
---
ChangeLog | 5 +
doc/examples_pkextract.dox | 20 ++--
doc/examples_pkmosaic.dox | 10 +-
src/algorithms/Filter2d.h | 2 +-
src/algorithms/StatFactory.h | 49 ++++-----
src/apps/Makefile.am | 3 +-
src/apps/pkcrop.cc | 9 +-
src/apps/pkdiff.cc | 2 +-
src/apps/pkextract.cc | 13 +--
src/apps/pkfilter.cc | 2 +-
src/apps/pkfilterdem.cc | 226 ++++++++++++++++++++++++++++++++++++++++
src/apps/pklas2img.cc | 238 +++++++++++++++++++++----------------------
src/apps/pkmosaic.cc | 123 +++++++++++++---------
13 files changed, 482 insertions(+), 220 deletions(-)
diff --git a/ChangeLog b/ChangeLog
index 4dd2a19..22b5b0f 100755
--- a/ChangeLog
+++ b/ChangeLog
@@ -239,3 +239,8 @@ version 2.5
renamed mask to nodata
- pkndvi
changed minmax option to min and max. Min and max values now apply to scaled values
+ - pklas2img
+ moved post filtering to pkfilterdem
+ changed some option names
+ - pkfilterdem
+ new utility (moved from pklas2ig)
diff --git a/doc/examples_pkextract.dox b/doc/examples_pkextract.dox
index 776e162..f9a0f51 100644
--- a/doc/examples_pkextract.dox
+++ b/doc/examples_pkextract.dox
@@ -2,36 +2,36 @@
\code
pkextract -i input.tif -s points.shp -o extracted.shp
\endcode
-extract all bands from input.tif to extracted.shp at pixel locations defined in points.shp.
+extract the points read in points.shp from input.tif. Create a new point vector file extracted.shp, where each point will contain an attribute for the individual input bands in input.tif.
\code
-pkextract -i input.tif -s points.shp -o extracted.shp -m mask.tif -msknodata 255
+pkextract -i input.tif -s points.shp -m mask.tif -msknodata 255 -o extracted.shp
\endcode
extract all bands from input.tif to extracted.shp at pixel locations defined in points.shp that have not a value 255 in mask.tif
\code
-pkextract -i input.tif -s locations.shp -o training.shp -c -1
+pkextract -i input.tif -s polygons.shp -o training.shp -r point
\endcode
-extract all classes (classes must have value different from 0) in input image input.tif to extracted.shp at pixel locations defined in points.shp
+extract all pixels from input.tif covered by the polygons in locations.shp. Each polygon can thus result in multiple point features with attributes for each input band. Write the extracted points to a point vector file training.shp.
\code
-pkextract -i input.tif -s polygons.shp -o extracted.shp -l -b 0
+pkextract -i input.tif -b 0 -s polygons.shp -r centroid -o extracted.shp -polygon
\endcode
-extract band 0 from input.tif to polygon vector file extracted.shp at centroids of polygons in polygons.shp.
+extract the first band from input.tif at the centroids of the polygons in vector filepolygons.shp. Assign the extracted point value to a new attribute of the polygon and write to the vector file extracted.shp.
\code
-pkextract -i input.tif -s polygons.shp -o extracted.shp -l -r 1 -b 1
+pkextract -i input.tif -b 1 -s polygons.shp -r mean -o extracted.shp -polygon
\endcode
-extract band 1 from input.tif to polygon vector file extracted.shp as means of polygons in polygons.shp.
+extract the mean values for the second band in input.tif covered by each polygon in polygons.shp. The mean values are written to a copy of the polygons in output vector file extracted.shp
\code
pkextract -i input.tif -s sample.tif -o extracted.shp -t 10 -c 1 -c 2 -c 3
\endcode
-extract all bands for a random sample of 10 percent of the pixels in sample.tif where sample.tif equals 1,2 or 3 (class values) to extracted.shp.
+Typical use where pixels are extracted based on a land cover map (sample.tif). Extract all bands for a random sample of 10 percent of the pixels in the land cover map sample.tif where the land cover classes are either 1,2 or 3 (class values). Write output to the point vector file extracted.shp.
\code
pkextract -i input.tif -s sample.tif -o extracted.shp -t -5000 -c 1
\endcode
-extract all bands for the first 5000 pixels (if available) where sample.tif equals 1 to extracted.shp.
+extract all bands for the first 5000 pixels encountered in sample.tif where pixels have a value equal to 1. Write output to point vector file extracted.shp.
diff --git a/doc/examples_pkmosaic.dox b/doc/examples_pkmosaic.dox
index 351ea9e..4be5f2c 100644
--- a/doc/examples_pkmosaic.dox
+++ b/doc/examples_pkmosaic.dox
@@ -5,14 +5,14 @@ pkmosaic -i input1.tif -i input2.tif -o output.tif
create mosaic from two input images. If images overlap, keep only last image (default rule)
\code
-pkmosaic -i input1.tif -i input2.tif -o output.tif --validBand 1 -srcnodata 255 -dstnodata 0
+pkmosaic -i input1.tif -i input2.tif -srcnodata 255 -bndnodata 1 -dstnodata 0 -o output.tif
\endcode
-create mosaic from two input images. Values of 255 in band 1 (starting from 0) are masked as invalid. Typically used when multi-band image contains cloud mask
+create mosaic from two input images. Values of 255 in band 1 (starting from 0) are masked as invalid. Typically used when second band of input image is a cloud mask
\code
-pkmosaic -i input1.tif -i input2.tif -o output.tif -cr maxndvi -rb 0 -rb 1 -srcnodata 255 -dstnodata 0
+pkmosaic -i input1.tif -i input2.tif -cr maxndvi -rb 0 -rb 1 -srcnodata 255 -bndnodata 0 -dstnodata 0 -o output.tif
\endcode
-create maximum NDVI (normalized difference vegetation index) composit. Values of 255 in band 0 (default) are masked as invalid and flagged as 0 if no other valid coverage. Typically used for (e.g., MODIS) images where red and near infrared spectral bands are stored in bands 0 and 1 respectively
+create maximum NDVI (normalized difference vegetation index) composit. Values of 255 in band 0 are masked as invalid and flagged as 0 if no other valid coverage. Typically used for (e.g., MODIS) images where red and near infrared spectral bands are stored in bands 0 and 1 respectively. In this particular case, a value of 255 in the first input band indicates a nodata value (e.g., cloud mask is coded within the data values).
\code
pkmosaic -i input1.tif -i input2.tif -i input3.tif -o output.tif -cr mean -w 0.75 -w 1.5 -w 0.75
@@ -22,5 +22,5 @@ create composite image using weighted mean: output=(3/4*input1+6/4*input2+3/4*in
\code
pkmosaic -i large.tif $(for IMAGE in *.tif;do pkinfo -i $IMAGE --cover $(pkinfo -i coverage.tif -bb);done) -cr median -min 0 -o output.tif
\endcode
-create median composit of all GTiff images found in current directory that cover (at least part of) coverage.tif. Values smaller or equal to 0 are set as nodatas 0 (default value for -dstnodata)
+create median composit of all GTiff images found in current directory that cover (at least part of) the image coverage.tif. Values smaller or equal to 0 are set as nodata 0 (default value for -dstnodata)
diff --git a/src/algorithms/Filter2d.h b/src/algorithms/Filter2d.h
index 272cb49..70f33cb 100644
--- a/src/algorithms/Filter2d.h
+++ b/src/algorithms/Filter2d.h
@@ -658,7 +658,7 @@ template<class T> unsigned long int Filter2d::morphology(const Vector2d<T>& inpu
break;
}
}
- output[y][x]=currentValue;//introduced due to hThrehold
+ output[y][x]=currentValue;//introduced due to hThreshold
if(currentMasked){
output[y][x]=currentValue;
}
diff --git a/src/algorithms/StatFactory.h b/src/algorithms/StatFactory.h
index c110c66..0b95b76 100644
--- a/src/algorithms/StatFactory.h
+++ b/src/algorithms/StatFactory.h
@@ -104,6 +104,7 @@ public:
gsl_rng_set(r,theSeed);
return r;
}
+ void getNodataValues(std::vector<double>& nodatav) const{nodatav=m_noDataValues;};
bool isNoData(double value) const{
if(m_noDataValues.empty())
return false;
@@ -119,7 +120,7 @@ public:
m_noDataValues=vnodata;
return m_noDataValues.size();
};
- static double getRandomValue(const gsl_rng* r, const std::string type, double a=0, double b=0){
+ double getRandomValue(const gsl_rng* r, const std::string type, double a=0, double b=0) const{
std::map<std::string, DISTRIBUTION_TYPE> m_distMap;
initDist(m_distMap);
double randValue=0;
@@ -161,26 +162,26 @@ public:
template<class T> double kurtosis(const std::vector<T>& v) const;
template<class T> void meanVar(const std::vector<T>& v, double& m1, double& v1) const;
template<class T1, class T2> void scale2byte(const std::vector<T1>& input, std::vector<T2>& output, unsigned char lbound=0, unsigned char ubound=255) const;
- template<class T> void 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=0.0, T &maximum=0.0, double sigma=0, const std::string &filename="");
- template<class T> void distribution(const std::vector<T>& input, std::vector<double>& output, int nbin, double sigma=0, const std::string &filename=""){distribution(input,input.begin(),input.end(),output,nbin,0,0,sigma,filename);};
- template<class T> void distribution2d(const std::vector<T>& inputX, const std::vector<T>& inputY, std::vector< std::vector<double> >& output, int nbin, T& minX=0, T& maxX=0, T& minY=0, T& maxY=0, double sigma=0, const std::string& filename="");
- template<class T> void 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);
- template<class T> void 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=10, T &minimum=0.0, T &maximum=0.0, const std::string &filename="");
- template<class T> void signature(const std::vector<T>& input, double& k, double& alpha, double& beta, double e);
- void signature(double m1, double m2, double& k, double& alpha, double& beta, double e);
- template<class T> void normalize(const std::vector<T>& input, std::vector<double>& output);
- template<class T> void normalize_pct(std::vector<T>& input);
+ template<class T> void 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=0.0, T &maximum=0.0, double sigma=0, const std::string &filename="") const;
+ template<class T> void distribution(const std::vector<T>& input, std::vector<double>& output, int nbin, double sigma=0, const std::string &filename="") const{distribution(input,input.begin(),input.end(),output,nbin,0,0,sigma,filename);};
+ template<class T> void distribution2d(const std::vector<T>& inputX, const std::vector<T>& inputY, std::vector< std::vector<double> >& output, int nbin, T& minX=0, T& maxX=0, T& minY=0, T& maxY=0, double sigma=0, const std::string& filename="") const;
+ template<class T> void 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) const;
+ template<class T> void 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=10, T &minimum=0.0, T &maximum=0.0, const std::string &filename="") const;
+ template<class T> void signature(const std::vector<T>& input, double& k, double& alpha, double& beta, double e) const;
+ void signature(double m1, double m2, double& k, double& alpha, double& beta, double e) const;
+ template<class T> void normalize(const std::vector<T>& input, std::vector<double>& output) const;
+ template<class T> void normalize_pct(std::vector<T>& input) const;
template<class T> double rmse(const std::vector<T>& x, const std::vector<T>& y) const;
template<class T> double correlation(const std::vector<T>& x, const std::vector<T>& y, int delay=0) const;
template<class T> double cross_correlation(const std::vector<T>& x, const std::vector<T>& y, int maxdelay, std::vector<T>& z) const;
template<class T> double linear_regression(const std::vector<T>& x, const std::vector<T>& y, double &c0, double &c1) const;
- template<class T> void 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=false);
- template<class T> void interpolateUp(const std::vector<double>& wavelengthIn, const std::vector< std::vector<T> >& input, const std::vector<double>& wavelengthOut, const std::string& type, std::vector< std::vector<T> >& output, bool verbose=false);
+ template<class T> void 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=false) const;
+ template<class T> void interpolateUp(const std::vector<double>& wavelengthIn, const std::vector< std::vector<T> >& input, const std::vector<double>& wavelengthOut, const std::string& type, std::vector< std::vector<T> >& output, bool verbose=false) const;
// template<class T> void interpolateUp(const std::vector< std::vector<T> >& input, std::vector< std::vector<T> >& output, double start, double end, double step, const gsl_interp_type* type);
// template<class T> void interpolateUp(const std::vector< std::vector<T> >& input, const std::vector<double>& wavelengthIn, std::vector< std::vector<T> >& output, std::vector<double>& wavelengthOut, double start, double end, double step, const gsl_interp_type* type);
- template<class T> void interpolateUp(const std::vector<T>& input, std::vector<T>& output, int nbin);
+ template<class T> void interpolateUp(const std::vector<T>& input, std::vector<T>& output, int nbin) const;
template<class T> void interpolateUp(double* input, int dim, std::vector<T>& output, int nbin);
- template<class T> void interpolateDown(const std::vector<T>& input, std::vector<T>& output, int nbin);
+ template<class T> void interpolateDown(const std::vector<T>& input, std::vector<T>& output, int nbin) const;
template<class T> void interpolateDown(double* input, int dim, std::vector<T>& output, int nbin);
private:
@@ -426,7 +427,7 @@ template<class T> inline double StatFactory::mean(const std::vector<T>& v) const
template<class T> inline T StatFactory::eraseNoData(std::vector<T>& v) const
{
- typename std::vector<T>::iterator it;
+ typename std::vector<T>::iterator it=v.begin();
while(it!=v.end()){
if(isNoData(*it))
v.erase(it);
@@ -568,7 +569,7 @@ template<class T1, class T2> void StatFactory::scale2byte(const std::vector<T1>&
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)
+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
{
double minValue=0;
double maxValue=0;
@@ -640,7 +641,7 @@ template<class T> void StatFactory::distribution(const std::vector<T>& input, t
}
}
-template<class T> void StatFactory::distribution2d(const std::vector<T>& inputX, const std::vector<T>& inputY, std::vector< std::vector<double> >& output, int nbin, T& minX, T& maxX, T& minY, T& maxY, double sigma, const std::string& filename)
+template<class T> void StatFactory::distribution2d(const std::vector<T>& inputX, const std::vector<T>& inputY, std::vector< std::vector<double> >& output, int nbin, T& minX, T& maxX, T& minY, T& maxY, double sigma, const std::string& filename) const
{
assert(inputX.size());
assert(inputX.size()==inputY.size());
@@ -720,7 +721,7 @@ template<class T> void StatFactory::distribution2d(const std::vector<T>& inputX
}
}
-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)
+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);
@@ -788,14 +789,14 @@ template<class T> void StatFactory::percentiles (const std::vector<T>& input, t
// }
// }
-template<class T> void StatFactory::signature(const std::vector<T>& input, double&k, double& alpha, double& beta, double e)
+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);
double m2=moment(input,2);
signature(m1,m2,k,alpha,beta,e);
}
-template<class T> void StatFactory::normalize(const std::vector<T>& input, std::vector<double>& output){
+template<class T> void StatFactory::normalize(const std::vector<T>& input, std::vector<double>& output) const{
double total=sum(input);
if(total){
output.resize(input.size());
@@ -806,7 +807,7 @@ template<class T> void StatFactory::normalize(const std::vector<T>& input, std::
output=input;
}
-template<class T> void StatFactory::normalize_pct(std::vector<T>& input){
+template<class T> void StatFactory::normalize_pct(std::vector<T>& input) const{
double total=sum(input);
if(total){
typename std::vector<T>::iterator it;
@@ -879,7 +880,7 @@ template<class T> double StatFactory::cross_correlation(const std::vector<T>& x,
//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)
-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){
+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());
assert(wavelengthOut.size());
@@ -950,7 +951,7 @@ template<class T> void StatFactory::interpolateUp(const std::vector<double>& wav
// gsl_interp_accel_free(acc);
// }
-template<class T> void StatFactory::interpolateUp(const std::vector<T>& input, std::vector<T>& output, int nbin)
+template<class T> void StatFactory::interpolateUp(const std::vector<T>& input, std::vector<T>& output, int nbin) const
{
assert(input.size());
assert(nbin);
@@ -990,7 +991,7 @@ template<class T> void StatFactory::interpolateUp(double* input, int dim, std::v
}
}
-template<class T> void StatFactory::interpolateDown(const std::vector<T>& input, std::vector<T>& output, int nbin)
+template<class T> void StatFactory::interpolateDown(const std::vector<T>& input, std::vector<T>& output, int nbin) const
{
assert(input.size());
assert(nbin);
diff --git a/src/apps/Makefile.am b/src/apps/Makefile.am
index 16e7ed1..2fff36f 100644
--- a/src/apps/Makefile.am
+++ b/src/apps/Makefile.am
@@ -6,7 +6,7 @@ LDADD = $(GSL_LIBS) $(GDAL_LDFLAGS) $(top_builddir)/src/algorithms/libalgorithms
###############################################################################
# the program to build and install (the names of the final binaries)
-bin_PROGRAMS = pkinfo pkcrop pkreclass pkgetmask pksetmask pkcreatect pkdumpimg pkdumpogr pksieve pkstatascii pkstatogr pkegcs pkextract pkfillnodata pkfilter pkenhance pkfilterascii pkdsm2shadow pkmosaic pkndvi pkpolygonize pkascii2img pkdiff pkclassify_svm pkfs_svm pkascii2ogr pkeditogr
+bin_PROGRAMS = pkinfo pkcrop pkreclass pkgetmask pksetmask pkcreatect pkdumpimg pkdumpogr pksieve pkstatascii pkstatogr pkegcs pkextract pkfillnodata pkfilter pkfilterdem pkenhance pkfilterascii pkdsm2shadow pkmosaic pkndvi pkpolygonize pkascii2img pkdiff pkclassify_svm pkfs_svm pkascii2ogr pkeditogr
# the program to build but not install (the names of the final binaries)
#noinst_PROGRAMS = pkxcorimg pkgeom
@@ -53,6 +53,7 @@ pkextract_SOURCES = pkextract.cc
pkfillnodata_SOURCES = pkfillnodata.cc
pkfilter_SOURCES = pkfilter.cc
pkfilter_LDADD = $(GSL_LIBS) $(AM_LDFLAGS) -lgsl -lgdal
+pkfilterdem_SOURCES = pkfilterdem.cc
pkenhance_SOURCES = pkenhance.cc
pkenhance_LDADD = $(AM_LDFLAGS) -lgdal
pkfilterascii_SOURCES = pkfilterascii.cc
diff --git a/src/apps/pkcrop.cc b/src/apps/pkcrop.cc
index 45d2a0e..899c415 100644
--- a/src/apps/pkcrop.cc
+++ b/src/apps/pkcrop.cc
@@ -137,8 +137,13 @@ int main(int argc, char *argv[])
ImgWriterGdal imgWriter;
//open input images to extract number of bands and spatial resolution
int ncropband=0;//total number of bands to write
- double dx=(dx_opt.size())? dx_opt[0]:0;
- double dy=(dy_opt.size())? dy_opt[0]:0;
+ double dx=0;
+ double dy=0;
+ if(dx_opt.size())
+ dx=dx_opt[0];
+ if(dy_opt.size())
+ dy=dy_opt[0];
+
for(int iimg=0;iimg<input_opt.size();++iimg){
imgReader.open(input_opt[iimg]);
if(dx_opt.empty()){
diff --git a/src/apps/pkdiff.cc b/src/apps/pkdiff.cc
index 448cc0d..d299fe5 100644
--- a/src/apps/pkdiff.cc
+++ b/src/apps/pkdiff.cc
@@ -45,7 +45,7 @@ int main(int argc, char *argv[])
Optionpk<string> labelref_opt("lr", "lref", "name of the reference label in case reference is shape file(default is label)", "label");
Optionpk<string> labelclass_opt("lc", "lclass", "name of the classified label in case output is shape file (default is class)", "class");
Optionpk<short> boundary_opt("bnd", "boundary", "boundary for selecting the sample (default: 1)", 1,1);
- Optionpk<bool> disc_opt("\0", "circular", "use circular disc kernel boundary)", false,1);
+ Optionpk<bool> disc_opt("circ", "circular", "use circular disc kernel boundary)", false,1);
Optionpk<bool> homogeneous_opt("hom", "homogeneous", "only take homogeneous regions into account", false,1);
Optionpk<string> option_opt("co", "co", "Creation option for output file. Multiple options can be specified.");
Optionpk<string> classname_opt("c", "class", "list of class names.");
diff --git a/src/apps/pkextract.cc b/src/apps/pkextract.cc
index ed911e7..eaac497 100644
--- a/src/apps/pkextract.cc
+++ b/src/apps/pkextract.cc
@@ -45,7 +45,7 @@ int main(int argc, char *argv[])
Optionpk<string> sample_opt("s", "sample", "Input sample vector file or class file (e.g. Corine CLC) if class option is set");
Optionpk<string> mask_opt("m", "mask", "Mask image file");
Optionpk<int> msknodata_opt("msknodata", "msknodata", "Mask value where image is invalid. If a single mask is used, more nodata values can be set. If more masks are used, use one value for each mask.", 1);
- Optionpk<int> class_opt("c", "class", "Class(es) to extract from input sample image. Use -1 to process every class in sample image, or leave empty to extract all valid data pixels from sample file");
+ Optionpk<int> class_opt("c", "class", "Class(es) to extract from input sample image. Leave empty to extract all valid data pixels from sample file");
Optionpk<string> output_opt("o", "output", "Output sample file (image file)");
Optionpk<string> ogrformat_opt("f", "f", "Output sample file format","ESRI Shapefile");
Optionpk<string> test_opt("test", "test", "Test sample file (use this option in combination with threshold<100 to create a training (output) and test set");
@@ -254,7 +254,7 @@ int main(int argc, char *argv[])
if(sampleIsRaster){
if(class_opt.empty()){
- std::cout << "Warning: no classes selected, if classes must be extracted, set to -1 for all classes using option -c -1" << std::endl;
+ // std::cout << "Warning: no classes selected, if a classes must be extracted, set to -1 for all classes using option -c -1" << std::endl;
ImgReaderGdal classReader;
ImgWriterOgr ogrWriter;
// if(verbose_opt[0]>1){
@@ -494,15 +494,10 @@ int main(int argc, char *argv[])
}
classReader.close();
nsample=writeBuffer.size();
- if(verbose_opt[0]){
+ if(verbose_opt[0])
std::cout << "total number of samples written: " << nsample << std::endl;
- if(nvalid.size()==class_opt.size()){
- for(int iclass=0;iclass<class_opt.size();++iclass)
- std::cout << "class " << class_opt[iclass] << " has " << nvalid[iclass] << " samples" << std::endl;
- }
- }
}
- else{//class_opt[0]!=0
+ else{//class_opt.size()!=0
assert(class_opt[0]);
// if(class_opt[0]){
assert(threshold_opt.size()==1||threshold_opt.size()==class_opt.size());
diff --git a/src/apps/pkfilter.cc b/src/apps/pkfilter.cc
index 0a0ce6c..0efcd31 100644
--- a/src/apps/pkfilter.cc
+++ b/src/apps/pkfilter.cc
@@ -40,7 +40,7 @@ int main(int argc,char **argv) {
Optionpk<std::string> input_opt("i","input","input image file");
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("c", "circular", "circular disc kernel for dilation and erosion", false);
+ Optionpk<bool> disc_opt("circ", "circular", "circular disc kernel for dilation and erosion", false);
// Optionpk<double> angle_opt("a", "angle", "angle used for directional filtering in dilation (North=0, East=90, South=180, West=270).");
Optionpk<std::string> method_opt("f", "filter", "filter function (median, var, min, max, sum, mean, dilate, erode, close, open, homog (central pixel must be identical to all other pixels within window), heterog, sobelx (horizontal edge detection), sobely (vertical edge detection), sobelxy (diagonal edge detection NE-SW),sobelyx (diagonal edge detection NW-SE), smooth, density, majority voting (only for classes), smoothnodata (smooth nodata values only) values, threshold local filtering [...]
Optionpk<std::string> resample_opt("r", "resampling-method", "Resampling method for shifting operation (near: nearest neighbour, bilinear: bi-linear interpolation).", "near");
diff --git a/src/apps/pkfilterdem.cc b/src/apps/pkfilterdem.cc
new file mode 100644
index 0000000..77f857e
--- /dev/null
+++ b/src/apps/pkfilterdem.cc
@@ -0,0 +1,226 @@
+/**********************************************************************
+pkfilterdem.cc: program to post filter raster images created with pklas2img
+Copyright (C) 2008-2014 Pieter Kempeneers
+
+This file is part of pktools
+
+pktools is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+pktools is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with pktools. If not, see <http://www.gnu.org/licenses/>.
+***********************************************************************/
+#include <assert.h>
+#include <iostream>
+#include <string>
+#include "base/Optionpk.h"
+#include "base/Vector2d.h"
+#include "algorithms/Filter2d.h"
+#include "imageclasses/ImgReaderGdal.h"
+#include "imageclasses/ImgWriterGdal.h"
+
+using namespace std;
+/*------------------
+ Main procedure
+ ----------------*/
+int main(int argc,char **argv) {
+ Optionpk<std::string> input_opt("i","input","input image file");
+ 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<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);
+ Optionpk<short> minChange_opt("minchange", "minchange", "Stop iterations when no more pixels are changed than this threshold.", 0);
+ Optionpk<std::string> otype_opt("ot", "otype", "Data type for output image ({Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/CInt16/CInt32/CFloat32/CFloat64}). Empty string: inherit type from input image","");
+ 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> verbose_opt("v", "verbose", "verbose mode if > 0", 0);
+
+ bool doProcess;//stop process when program was invoked with help option (-h --help)
+ try{
+ doProcess=input_opt.retrieveOption(argc,argv);
+ output_opt.retrieveOption(argc,argv);
+ tmpdir_opt.retrieveOption(argc,argv);
+ disc_opt.retrieveOption(argc,argv);
+ postFilter_opt.retrieveOption(argc,argv);
+ dim_opt.retrieveOption(argc,argv);
+ maxSlope_opt.retrieveOption(argc,argv);
+ hThreshold_opt.retrieveOption(argc,argv);
+ minChange_opt.retrieveOption(argc,argv);
+ otype_opt.retrieveOption(argc,argv);
+ oformat_opt.retrieveOption(argc,argv);
+ colorTable_opt.retrieveOption(argc,argv);
+ verbose_opt.retrieveOption(argc,argv);
+ }
+ catch(string predefinedString){
+ std::cout << predefinedString << std::endl;
+ exit(0);
+ }
+ if(!doProcess){
+ std::cout << "short option -h shows basic options only, use long option --help to show all options" << std::endl;
+ exit(0);//help was invoked, stop processing
+ }
+
+ ImgReaderGdal input;
+ ImgWriterGdal outputWriter;
+ if(input_opt.empty()){
+ cerr << "Error: no input file selected, use option -i" << endl;
+ exit(1);
+ }
+ if(output_opt.empty()){
+ cerr << "Error: no outputWriter file selected, use option -o" << endl;
+ exit(1);
+ }
+ if(postFilter_opt.empty()){
+ cerr << "Error: no filter selected, use option -f" << endl;
+ exit(1);
+ }
+ input.open(input_opt[0]);
+ GDALDataType theType=GDT_Unknown;
+ if(verbose_opt[0])
+ cout << "possible output data types: ";
+ for(int iType = 0; iType < GDT_TypeCount; ++iType){
+ if(verbose_opt[0])
+ cout << " " << GDALGetDataTypeName((GDALDataType)iType);
+ if( GDALGetDataTypeName((GDALDataType)iType) != NULL
+ && EQUAL(GDALGetDataTypeName((GDALDataType)iType),
+ otype_opt[0].c_str()))
+ theType=(GDALDataType) iType;
+ }
+ if(theType==GDT_Unknown)
+ theType=input.getDataType();
+
+ if(verbose_opt[0])
+ std::cout << std::endl << "Output pixel type: " << GDALGetDataTypeName(theType) << endl;
+
+ string imageType=input.getImageType();
+ if(oformat_opt.size())
+ imageType=oformat_opt[0];
+
+ if(option_opt.findSubstring("INTERLEAVE=")==option_opt.end()){
+ string theInterleave="INTERLEAVE=";
+ theInterleave+=input.getInterleave();
+ option_opt.push_back(theInterleave);
+ }
+
+ if(verbose_opt[0])
+ cout << "opening output file " << output_opt[0] << endl;
+ outputWriter.open(output_opt[0],input.nrOfCol(),input.nrOfRow(),1,theType,imageType,option_opt);
+ //set projection
+ outputWriter.setProjection(input.getProjection());
+ outputWriter.copyGeoTransform(input);
+ if(colorTable_opt.size())
+ outputWriter.setColorTable(colorTable_opt[0]);
+
+ Vector2d<double> inputData(input.nrOfRow(),input.nrOfCol());
+ Vector2d<double> outputData(outputWriter.nrOfRow(),outputWriter.nrOfCol());
+ Vector2d<double> tmpData(outputWriter.nrOfRow(),outputWriter.nrOfCol());
+ input.readDataBlock(inputData,GDT_Float64,0,inputData.nCols()-1,0,inputData.nRows()-1);
+
+ //apply post filter
+ std::cout << "Applying post processing filter: " << postFilter_opt[0] << std::endl;
+
+ // const char* pszMessage;
+ // void* pProgressArg=NULL;
+ // GDALProgressFunc pfnProgress=GDALTermProgress;
+ // double progress=0;
+ // pfnProgress(progress,pszMessage,pProgressArg);
+
+ //make sure dim_opt contains initial [0] and maximum [1] kernel sizes in this order
+ if(dim_opt.size()<2)
+ dim_opt.insert(dim_opt.begin(),3);
+ if(dim_opt[0]>dim_opt[1]){
+ dim_opt.insert(dim_opt.begin(),dim_opt[1]);
+ dim_opt.erase(dim_opt.begin()+2);
+ }
+ 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);
+ 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);
+ inputData=outputData;
+ dim+=2;//change from theory: originally double cellCize
+ std::cout << "iteration " << iteration << ": " << nchange << " pixels changed" << std::endl;
+ ++iteration;
+ }
+ }
+ else if(postFilter_opt[0]=="promorph"){//todo: instead of number of iterations, define a max dim size
+ //Progressive morphological filter tgrs2003_zhang vol41 pp 872-882
+ //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]){
+ std::cout << "iteration " << iteration << " with window size " << dim << " and dh_max: " << hThreshold << std::endl;
+ try{
+ nchange=theFilter.morphology(inputData,outputData,"erode",dim,dim,disc_opt[0],hThreshold);
+ theFilter.morphology(outputData,inputData,"dilate",dim,dim,disc_opt[0],hThreshold);
+ theFilter.doit(inputData,outputData,"median",dim,dim,1,disc_opt[0]);
+ inputData=outputData;
+ }
+ catch(std::string errorString){
+ cout << errorString << endl;
+ exit(1);
+ }
+ int newdim=(dim==1)? 3: 2*(dim-1)+1;
+ hThreshold=hThreshold_opt[0]+maxSlope_opt[0]*(newdim-dim)*input.getDeltaX();
+ dim=newdim;
+ if(hThreshold_opt.size()>1){
+ if(hThreshold>hThreshold_opt[1])
+ hThreshold=hThreshold_opt[1];
+ }
+ std::cout << "iteration " << iteration << ": " << nchange << " pixels changed" << std::endl;
+ ++iteration;
+ }
+ }
+ 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]);
+ outputData=inputData;
+ }
+ catch(std::string errorString){
+ cout << errorString << endl;
+ exit(1);
+ }
+ }
+ 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]);
+ }
+ catch(std::string errorString){
+ cout << errorString << endl;
+ exit(1);
+ }
+ }
+ //write outputData to outputWriter
+ outputWriter.writeDataBlock(outputData,GDT_Float64,0,outputData.nCols()-1,0,outputData.nRows()-1);
+
+ // progress=1;
+ // pfnProgress(progress,pszMessage,pProgressArg);
+ input.close();
+ outputWriter.close();
+ return 0;
+}
diff --git a/src/apps/pklas2img.cc b/src/apps/pklas2img.cc
index 22b6d8e..77aa6a0 100644
--- a/src/apps/pklas2img.cc
+++ b/src/apps/pklas2img.cc
@@ -35,21 +35,21 @@ int main(int argc,char **argv) {
Optionpk<bool> disc_opt("circ", "circular", "circular disc kernel for dilation and erosion", false);
Optionpk<double> maxSlope_opt("s", "maxSlope", "Maximum slope used for morphological filtering", 0.0);
Optionpk<double> hThreshold_opt("ht", "maxHeight", "initial and maximum height threshold for progressive morphological filtering (e.g., -ht 0.2 -ht 2.5)", 0.2);
- Optionpk<short> maxIter_opt("\0", "maxIter", "Maximum number of iterations in post filter", 100.0);
- Optionpk<short> nbin_opt("nb", "nbin", "Number of percentile bins for calculating profile (=number of output bands)", 10.0);
- Optionpk<unsigned short> returns_opt("r", "returns", "number(s) of returns to include");
- Optionpk<unsigned short> classes_opt("c", "classes", "classes to keep: 0 (created, never classified), 1 (unclassified), 2 (ground), 3 (low vegetation), 4 (medium vegetation), 5 (high vegetation), 6 (building), 7 (low point, noise), 8 (model key-point), 9 (water), 10 (reserved), 11 (reserved), 12 (overlap)");
+ Optionpk<short> maxIter_opt("maxit", "maxit", "Maximum number of iterations in post filter", 5);
+ Optionpk<short> nbin_opt("nbin", "nbin", "Number of percentile bins for calculating profile (=number of output bands)", 10.0);
+ Optionpk<unsigned short> returns_opt("ret", "ret", "number(s) of returns to include");
+ Optionpk<unsigned short> classes_opt("class", "class", "classes to keep: 0 (created, never classified), 1 (unclassified), 2 (ground), 3 (low vegetation), 4 (medium vegetation), 5 (high vegetation), 6 (building), 7 (low point, noise), 8 (model key-point), 9 (water), 10 (reserved), 11 (reserved), 12 (overlap)");
Optionpk<string> composite_opt("comp", "comp", "composite for multiple points in cell (min, max, median, mean, sum, first, last, profile, number (point density)). Last: overwrite cells with latest point", "last");
Optionpk<string> filter_opt("fir", "filter", "filter las points (last,single,multiple,all).", "all");
- Optionpk<string> postFilter_opt("pf", "pfilter", "post processing filter (etew_min,promorph (progressive morphological filter),bunting (adapted promorph),open,close,none).", "none");
- Optionpk<short> dimx_opt("\0", "dimX", "Dimension X of postFilter", 3);
- Optionpk<short> dimy_opt("\0", "dimY", "Dimension Y of postFilter", 3);
+ // Optionpk<string> postFilter_opt("pf", "pfilter", "post processing filter (etew_min,promorph (progressive morphological filter),bunting (adapted promorph),open,close,none).", "none");
+ // Optionpk<short> dimx_opt("dimx", "dimx", "Dimension X of postFilter", 3);
+ // Optionpk<short> dimy_opt("dimy", "dimy", "Dimension Y of postFilter", 3);
Optionpk<string> output_opt("o", "output", "Output image file");
Optionpk<string> projection_opt("a_srs", "a_srs", "assign the projection for the output file in epsg code, e.g., epsg:3035 for European LAEA projection");
- Optionpk<double> ulx_opt("\0", "ulx", "Upper left x value bounding box (in geocoordinates if georef is true). 0 is read from input file", 0.0);
- Optionpk<double> uly_opt("\0", "uly", "Upper left y value bounding box (in geocoordinates if georef is true). 0 is read from input file", 0.0);
- Optionpk<double> lrx_opt("\0", "lrx", "Lower right x value bounding box (in geocoordinates if georef is true). 0 is read from input file", 0.0);
- Optionpk<double> lry_opt("\0", "lry", "Lower right y value bounding box (in geocoordinates if georef is true). 0 is read from input file", 0.0);
+ Optionpk<double> ulx_opt("ulx", "ulx", "Upper left x value bounding box (in geocoordinates if georef is true). 0 is read from input file", 0.0);
+ Optionpk<double> uly_opt("uly", "uly", "Upper left y value bounding box (in geocoordinates if georef is true). 0 is read from input file", 0.0);
+ Optionpk<double> lrx_opt("lrx", "lrx", "Lower right x value bounding box (in geocoordinates if georef is true). 0 is read from input file", 0.0);
+ Optionpk<double> lry_opt("lry", "lry", "Lower right y value bounding box (in geocoordinates if georef is true). 0 is read from input file", 0.0);
Optionpk<string> otype_opt("ot", "otype", "Data type for output image ({Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/CInt16/CInt32/CFloat32/CFloat64}). Empty string: inherit type from input image", "Byte");
Optionpk<string> oformat_opt("of", "oformat", "Output image format (see also gdal_translate). Empty string: inherit from input image", "GTiff");
Optionpk<string> option_opt("co", "co", "Creation option for output file. Multiple options can be specified.");
@@ -74,9 +74,9 @@ int main(int argc,char **argv) {
classes_opt.retrieveOption(argc,argv);
composite_opt.retrieveOption(argc,argv);
filter_opt.retrieveOption(argc,argv);
- postFilter_opt.retrieveOption(argc,argv);
- dimx_opt.retrieveOption(argc,argv);
- dimy_opt.retrieveOption(argc,argv);
+ // postFilter_opt.retrieveOption(argc,argv);
+ // dimx_opt.retrieveOption(argc,argv);
+ // dimy_opt.retrieveOption(argc,argv);
output_opt.retrieveOption(argc,argv);
projection_opt.retrieveOption(argc,argv);
ulx_opt.retrieveOption(argc,argv);
@@ -190,7 +190,7 @@ int main(int argc,char **argv) {
std::cout << setprecision(12) << "--ulx=" << minULX << " --uly=" << maxULY << " --lrx=" << maxLRX << " --lry=" << minLRY << std::endl;
std::cout << "total number of points before filtering: " << totalPoints << std::endl;
std::cout << "filter set to " << filter_opt[0] << std::endl;
- std::cout << "postFilter set to " << postFilter_opt[0] << std::endl;
+ // std::cout << "postFilter set to " << postFilter_opt[0] << std::endl;
}
int ncol=ceil(maxLRX-minULX)/dx_opt[0];//number of columns in outputGrid
int nrow=ceil(maxULY-minLRY)/dy_opt[0];//number of rows in outputGrid
@@ -294,16 +294,14 @@ int main(int argc,char **argv) {
outputWriter.geo2image(thePoint.GetX(),thePoint.GetY(),dcol,drow);
int icol=static_cast<int>(dcol);
int irow=static_cast<int>(drow);
- // //test
- // irow+=1;
if(irow<0||irow>=nrow){
- //test
- cout << "Error: thePoint.GetX(),thePoint.GetY(),dcol,drow" << thePoint.GetX() << ", " << thePoint.GetY() << ", " << dcol << ", " << drow << endl;
+ // //test
+ // cout << "Error: thePoint.GetX(),thePoint.GetY(),dcol,drow" << thePoint.GetX() << ", " << thePoint.GetY() << ", " << dcol << ", " << drow << endl;
continue;
}
if(icol<0||icol>=ncol){
- //test
- cout << "Error: thePoint.GetX(),thePoint.GetY(),dcol,drow" << thePoint.GetX() << ", " << thePoint.GetY() << ", " << dcol << ", " << drow << endl;
+ // //test
+ // cout << "Error: thePoint.GetX(),thePoint.GetY(),dcol,drow" << thePoint.GetX() << ", " << thePoint.GetY() << ", " << dcol << ", " << drow << endl;
continue;
}
assert(irow>=0);
@@ -338,11 +336,11 @@ int main(int argc,char **argv) {
pfnProgress(progress,pszMessage,pProgressArg);
statfactory::StatFactory stat;
//fill inputData in outputData
- if(composite_opt[0]=="profile"){
- assert(postFilter_opt[0]=="none");
+ // if(composite_opt[0]=="profile"){
+ // assert(postFilter_opt[0]=="none");
// for(int iband=0;iband<nband;++iband)
// outputProfile[iband].resize(nrow,ncol);
- }
+ // }
for(int irow=0;irow<nrow;++irow){
if(composite_opt[0]=="number")
continue;//outputData already set
@@ -415,102 +413,102 @@ int main(int argc,char **argv) {
pfnProgress(progress,pszMessage,pProgressArg);
inputData.clear();//clean up memory
//apply post filter
- std::cout << "Applying post processing filter: " << postFilter_opt[0] << std::endl;
- if(postFilter_opt[0]=="etew_min"){
- if(composite_opt[0]!="min")
- std::cout << "Warning: composite option is not set to min!" << std::endl;
- //Elevation Threshold with Expand Window (ETEW) Filter (p.73 frmo Airborne LIDAR Data Processing and Analysis Tools ALDPAT 1.0)
- //first iteration is performed assuming only minima are selected using options -fir all -c min
- unsigned long int nchange=1;
- //increase cells and thresholds until no points from the previous iteration are discarded.
- int dimx=dimx_opt[0];
- int dimy=dimy_opt[0];
- filter2d::Filter2d morphFilter;
- // morphFilter.setNoValue(0);
- Vector2d<float> currentOutput=outputData;
- int iteration=1;
- while(nchange&&iteration<maxIter_opt[0]){
- double hThreshold=maxSlope_opt[0]*dimx;
- Vector2d<float> newOutput;
- nchange=morphFilter.morphology(currentOutput,newOutput,"erode",dimx,dimy,disc_opt[0],hThreshold);
- currentOutput=newOutput;
- dimx+=2;//change from theory: originally double cellCize
- dimy+=2;//change from theory: originally double cellCize
- std::cout << "iteration " << iteration << ": " << nchange << " pixels changed" << std::endl;
- ++iteration;
- }
- outputData=currentOutput;
- }
- else if(postFilter_opt[0]=="promorph"||postFilter_opt[0]=="bunting"){
- if(composite_opt[0]!="min")
- std::cout << "Warning: composite option is not set to min!" << std::endl;
- assert(hThreshold_opt.size()>1);
- //Progressive morphological filter tgrs2003_zhang vol41 pp 872-882
- //first iteration is performed assuming only minima are selected using options -fir all -c min
- //increase cells and thresholds until no points from the previous iteration are discarded.
- int dimx=dimx_opt[0];
- int dimy=dimy_opt[0];
- filter2d::Filter2d theFilter;
- // theFilter.setNoValue(0);
- Vector2d<float> currentOutput=outputData;
- double hThreshold=hThreshold_opt[0];
- int iteration=1;
- while(iteration<maxIter_opt[0]){
- std::cout << "iteration " << iteration << " with window size " << dimx << " and dh_max: " << hThreshold << std::endl;
- Vector2d<float> newOutput;
- try{
- theFilter.morphology(outputData,currentOutput,"erode",dimx,dimy,disc_opt[0],hThreshold);
- theFilter.morphology(currentOutput,outputData,"dilate",dimx,dimy,disc_opt[0],hThreshold);
- if(postFilter_opt[0]=="bunting"){
- theFilter.doit(outputData,currentOutput,"median",dimx,dimy,1,disc_opt[0]);
- outputData=currentOutput;
- }
- }
- catch(std::string errorString){
- cout << errorString << endl;
- exit(1);
- }
- int newdimx=(dimx==1)? 3: 2*(dimx-1)+1;
- int newdimy=(dimx==1)? 3: 2*(dimy-1)+1;//from PE&RS vol 71 pp313-324
- hThreshold=hThreshold_opt[0]+maxSlope_opt[0]*(newdimx-dimx)*dx_opt[0];
- dimx=newdimx;
- dimy=newdimy;
- if(hThreshold>hThreshold_opt[1])
- hThreshold=hThreshold_opt[1];
- ++iteration;
- }
- outputData=currentOutput;
- }
- else if(postFilter_opt[0]=="open"){
- if(composite_opt[0]!="min")
- std::cout << "Warning: composite option is not set to min!" << std::endl;
- filter2d::Filter2d morphFilter;
- // morphFilter.setNoValue(0);
- Vector2d<float> filterInput=outputData;
- try{
- morphFilter.morphology(outputData,filterInput,"erode",dimx_opt[0],dimy_opt[0],disc_opt[0],maxSlope_opt[0]);
- morphFilter.morphology(filterInput,outputData,"dilate",dimx_opt[0],dimy_opt[0],disc_opt[0],maxSlope_opt[0]);
- }
- catch(std::string errorString){
- cout << errorString << endl;
- exit(1);
- }
- }
- else if(postFilter_opt[0]=="close"){
- if(composite_opt[0]!="max")
- std::cout << "Warning: composite option is not set to max!" << std::endl;
- filter2d::Filter2d morphFilter;
- // morphFilter.setNoValue(0);
- Vector2d<float> filterInput=outputData;
- try{
- morphFilter.morphology(outputData,filterInput,"dilate",dimx_opt[0],dimy_opt[0],disc_opt[0],maxSlope_opt[0]);
- morphFilter.morphology(filterInput,outputData,"erode",dimx_opt[0],dimy_opt[0],disc_opt[0],maxSlope_opt[0]);
- }
- catch(std::string errorString){
- cout << errorString << endl;
- exit(1);
- }
- }
+ // std::cout << "Applying post processing filter: " << postFilter_opt[0] << std::endl;
+ // if(postFilter_opt[0]=="etew_min"){
+ // if(composite_opt[0]!="min")
+ // std::cout << "Warning: composite option is not set to min!" << std::endl;
+ // //Elevation Threshold with Expand Window (ETEW) Filter (p.73 frmo 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
+ // unsigned long int nchange=1;
+ // //increase cells and thresholds until no points from the previous iteration are discarded.
+ // int dimx=dimx_opt[0];
+ // int dimy=dimy_opt[0];
+ // filter2d::Filter2d morphFilter;
+ // // morphFilter.setNoValue(0);
+ // Vector2d<float> currentOutput=outputData;
+ // int iteration=1;
+ // while(nchange&&iteration<=maxIter_opt[0]){
+ // double hThreshold=maxSlope_opt[0]*dimx;
+ // Vector2d<float> newOutput;
+ // nchange=morphFilter.morphology(currentOutput,newOutput,"erode",dimx,dimy,disc_opt[0],hThreshold);
+ // currentOutput=newOutput;
+ // dimx+=2;//change from theory: originally double cellCize
+ // dimy+=2;//change from theory: originally double cellCize
+ // std::cout << "iteration " << iteration << ": " << nchange << " pixels changed" << std::endl;
+ // ++iteration;
+ // }
+ // outputData=currentOutput;
+ // }
+ // else if(postFilter_opt[0]=="promorph"||postFilter_opt[0]=="bunting"){
+ // if(composite_opt[0]!="min")
+ // std::cout << "Warning: composite option is not set to min!" << std::endl;
+ // assert(hThreshold_opt.size()>1);
+ // //Progressive morphological filter tgrs2003_zhang vol41 pp 872-882
+ // //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 dimx=dimx_opt[0];
+ // int dimy=dimy_opt[0];
+ // filter2d::Filter2d theFilter;
+ // // theFilter.setNoValue(0);
+ // Vector2d<float> currentOutput=outputData;
+ // double hThreshold=hThreshold_opt[0];
+ // int iteration=1;
+ // while(iteration<=maxIter_opt[0]){
+ // std::cout << "iteration " << iteration << " with window size " << dimx << " and dh_max: " << hThreshold << std::endl;
+ // Vector2d<float> newOutput;
+ // try{
+ // theFilter.morphology(outputData,currentOutput,"erode",dimx,dimy,disc_opt[0],hThreshold);
+ // theFilter.morphology(currentOutput,outputData,"dilate",dimx,dimy,disc_opt[0],hThreshold);
+ // if(postFilter_opt[0]=="bunting"){
+ // theFilter.doit(outputData,currentOutput,"median",dimx,dimy,1,disc_opt[0]);
+ // outputData=currentOutput;
+ // }
+ // }
+ // catch(std::string errorString){
+ // cout << errorString << endl;
+ // exit(1);
+ // }
+ // int newdimx=(dimx==1)? 3: 2*(dimx-1)+1;
+ // int newdimy=(dimx==1)? 3: 2*(dimy-1)+1;//from PE&RS vol 71 pp313-324
+ // hThreshold=hThreshold_opt[0]+maxSlope_opt[0]*(newdimx-dimx)*dx_opt[0];
+ // dimx=newdimx;
+ // dimy=newdimy;
+ // if(hThreshold>hThreshold_opt[1])
+ // hThreshold=hThreshold_opt[1];
+ // ++iteration;
+ // }
+ // outputData=currentOutput;
+ // }
+ // else if(postFilter_opt[0]=="open"){
+ // if(composite_opt[0]!="min")
+ // std::cout << "Warning: composite option is not set to min!" << std::endl;
+ // filter2d::Filter2d morphFilter;
+ // // morphFilter.setNoValue(0);
+ // Vector2d<float> filterInput=outputData;
+ // try{
+ // morphFilter.morphology(outputData,filterInput,"erode",dimx_opt[0],dimy_opt[0],disc_opt[0],maxSlope_opt[0]);
+ // morphFilter.morphology(filterInput,outputData,"dilate",dimx_opt[0],dimy_opt[0],disc_opt[0],maxSlope_opt[0]);
+ // }
+ // catch(std::string errorString){
+ // cout << errorString << endl;
+ // exit(1);
+ // }
+ // }
+ // else if(postFilter_opt[0]=="close"){
+ // if(composite_opt[0]!="max")
+ // std::cout << "Warning: composite option is not set to max!" << std::endl;
+ // filter2d::Filter2d morphFilter;
+ // // morphFilter.setNoValue(0);
+ // Vector2d<float> filterInput=outputData;
+ // try{
+ // morphFilter.morphology(outputData,filterInput,"dilate",dimx_opt[0],dimy_opt[0],disc_opt[0],maxSlope_opt[0]);
+ // morphFilter.morphology(filterInput,outputData,"erode",dimx_opt[0],dimy_opt[0],disc_opt[0],maxSlope_opt[0]);
+ // }
+ // catch(std::string errorString){
+ // cout << errorString << endl;
+ // exit(1);
+ // }
+ // }
if(composite_opt[0]!="profile"){
//write output file
std::cout << "writing output raster file" << std::endl;
diff --git a/src/apps/pkmosaic.cc b/src/apps/pkmosaic.cc
index 39fe661..4e3e4aa 100644
--- a/src/apps/pkmosaic.cc
+++ b/src/apps/pkmosaic.cc
@@ -38,28 +38,28 @@ int main(int argc, char *argv[])
{
Optionpk<string> input_opt("i", "input", "Input image file(s). If input contains multiple images, a multi-band output is created");
Optionpk<string> output_opt("o", "output", "Output image file");
- Optionpk<string> projection_opt("a_srs", "a_srs", "Override the projection for the output file (leave blank to copy from input file, use epsg:3035 to use European projection and force to European grid");
+ Optionpk<string> projection_opt("a_srs", "a_srs", "Override the spatial reference for the output file (leave blank to copy from input file, use epsg:3035 to use European projection and force to European grid");
Optionpk<string> extent_opt("e", "extent", "get boundary from extent from polygons in vector file");
Optionpk<double> ulx_opt("ulx", "ulx", "Upper left x value bounding box", 0.0);
Optionpk<double> uly_opt("uly", "uly", "Upper left y value bounding box", 0.0);
Optionpk<double> lrx_opt("lrx", "lrx", "Lower right x value bounding box", 0.0);
Optionpk<double> lry_opt("lry", "lry", "Lower right y value bounding box", 0.0);
- Optionpk<double> dx_opt("dx", "dx", "Output resolution in x (in meter) (0.0: keep original resolution)", 0.0);
- Optionpk<double> dy_opt("dy", "dy", "Output resolution in y (in meter) (0.0: keep original resolution)", 0.0);
+ Optionpk<double> dx_opt("dx", "dx", "Output resolution in x (in meter) (empty: keep original resolution)");
+ Optionpk<double> dy_opt("dy", "dy", "Output resolution in y (in meter) (empty: keep original resolution)");
Optionpk<int> band_opt("b", "band", "band index(es) to crop (leave empty if all bands must be retained)");
Optionpk<string> otype_opt("ot", "otype", "Data type for output image ({Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/CInt16/CInt32/CFloat32/CFloat64}). Empty string: inherit type from input image", "");
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)");
Optionpk<string> option_opt("co", "co", "Creation option for output file. Multiple options can be specified.");
- Optionpk<double> dstnodata_opt("dstnodata", "dstnodata", "nodata value to put in image if out of bounds.", 0);
Optionpk<unsigned short> resample_opt("r", "resample", "Resampling method (0: nearest neighbour, 1: bi-linear interpolation).", 0);
Optionpk<string> description_opt("d", "description", "Set image description");
Optionpk<string> crule_opt("cr", "crule", "Composite rule for mosaic (overwrite, maxndvi, maxband, minband, mean, mode (only for byte images), median, sum", "overwrite");
Optionpk<int> ruleBand_opt("rb", "rband", "band index used for the rule (for ndvi, use --ruleBand=redBand --ruleBand=nirBand", 0);
- Optionpk<int> validBand_opt("vb", "validBand", "valid band index(es)", 0);
- Optionpk<double> srcnodata_opt("srcnodata", "srcnodata", "invalid value for valid band", 0);
- Optionpk<double> minValue_opt("min", "min", "flag values smaller or equal to this value as invalid.", -99999999);
- Optionpk<double> maxValue_opt("max", "max", "flag values larger or equal to this value as invalid.", 99999999);
+ Optionpk<double> srcnodata_opt("srcnodata", "srcnodata", "invalid value for input image", 0);
+ Optionpk<int> bndnodata_opt("bndnodata", "bndnodata", "Bands in input image to check if pixel is valid (used for srcnodata, min and max options)", 0);
+ Optionpk<double> dstnodata_opt("dstnodata", "dstnodata", "nodata value to put in output image if not valid or out of bounds.", 0);
+ Optionpk<double> minValue_opt("min", "min", "flag values smaller or equal to this value as invalid.");
+ Optionpk<double> maxValue_opt("max", "max", "flag values larger or equal to this value as invalid.");
Optionpk<bool> file_opt("file", "file", "write number of observations for each pixels as additional layer in mosaic", false);
Optionpk<short> weight_opt("w", "weight", "Weights (type: short) for the mosaic, use one weight for each input file in same order as input files are provided). Use value 1 for equal weights.", 1);
Optionpk<short> class_opt("c", "class", "classes for multi-band output image: each band represents the number of observations for one specific class. Use value 0 for no multi-band output image).", 0);
@@ -87,7 +87,7 @@ int main(int argc, char *argv[])
description_opt.retrieveOption(argc,argv);
crule_opt.retrieveOption(argc,argv);
ruleBand_opt.retrieveOption(argc,argv);
- validBand_opt.retrieveOption(argc,argv);
+ bndnodata_opt.retrieveOption(argc,argv);
srcnodata_opt.retrieveOption(argc,argv);
minValue_opt.retrieveOption(argc,argv);
maxValue_opt.retrieveOption(argc,argv);
@@ -119,18 +119,22 @@ int main(int argc, char *argv[])
cruleMap["median"]=crule::median;
cruleMap["sum"]=crule::sum;
- while(srcnodata_opt.size()<validBand_opt.size())
+ while(srcnodata_opt.size()<bndnodata_opt.size())
srcnodata_opt.push_back(srcnodata_opt[0]);
- while(validBand_opt.size()<srcnodata_opt.size())
- validBand_opt.push_back(validBand_opt[0]);
- while(minValue_opt.size()<validBand_opt.size())
- minValue_opt.push_back(minValue_opt[0]);
- while(validBand_opt.size()<minValue_opt.size())
- validBand_opt.push_back(validBand_opt[0]);
- while(maxValue_opt.size()<validBand_opt.size())
- maxValue_opt.push_back(maxValue_opt[0]);
- while(validBand_opt.size()<maxValue_opt.size())
- validBand_opt.push_back(validBand_opt[0]);
+ while(bndnodata_opt.size()<srcnodata_opt.size())
+ bndnodata_opt.push_back(bndnodata_opt[0]);
+ if(minValue_opt.size()){
+ while(minValue_opt.size()<bndnodata_opt.size())
+ minValue_opt.push_back(minValue_opt[0]);
+ while(bndnodata_opt.size()<minValue_opt.size())
+ bndnodata_opt.push_back(bndnodata_opt[0]);
+ }
+ if(maxValue_opt.size()){
+ while(maxValue_opt.size()<bndnodata_opt.size())
+ maxValue_opt.push_back(maxValue_opt[0]);
+ while(bndnodata_opt.size()<maxValue_opt.size())
+ bndnodata_opt.push_back(bndnodata_opt[0]);
+ }
RESAMPLE theResample;
switch(resample_opt[0]){
case(BILINEAR):
@@ -176,8 +180,8 @@ int main(int argc, char *argv[])
cout << "Output pixel type: " << GDALGetDataTypeName(theType) << endl;
}
- double dx=dx_opt[0];
- double dy=dy_opt[0];
+ double dx=0;
+ double dy=0;
//get bounding box from extentReader if defined
ImgReaderOgr extentReader;
if(extent_opt.size()){
@@ -274,13 +278,8 @@ int main(int argc, char *argv[])
for(int iband=0;iband<nband;++iband)
bands[iband]=iband;
}
- assert(validBand_opt.size()==minValue_opt.size());
- assert(validBand_opt.size()==maxValue_opt.size());
- for(int iband=0;iband<validBand_opt.size();++iband){
- assert(validBand_opt[iband]>=0&&validBand_opt[iband]<nband);
- if(verbose_opt[0]){
- cout << "band " << validBand_opt[iband] << " is valid in ] " << minValue_opt[iband] << " , " << maxValue_opt[iband] << " [" << endl;
- }
+ for(int iband=0;iband<bndnodata_opt.size();++iband){
+ assert(bndnodata_opt[iband]>=0&&bndnodata_opt[iband]<nband);
}
//if output type not set, get type from input image
if(theType==GDT_Unknown){
@@ -304,10 +303,14 @@ int main(int argc, char *argv[])
maxULY=theULY;
minULX=theULX;
minLRY=theLRY;
- if(!dx||!dy){
+ if(dx_opt.size())
+ dx=dx_opt[0];
+ else
dx=imgReader.getDeltaX();
+ if(dy_opt.size())
+ dy=dy_opt[0];
+ else
dy=imgReader.getDeltaY();
- }
// imgReader.getMagicPixel(magic_x,magic_y);
init=true;
}
@@ -543,25 +546,53 @@ int main(int argc, char *argv[])
lowerCol=0;
if(upperCol>=imgReader.nrOfCol())
upperCol=imgReader.nrOfCol()-1;
- for(int vband=0;vband<validBand_opt.size();++vband){
- val_new=(readCol-0.5-lowerCol)*readBuffer[validBand_opt[vband]][upperCol-startCol]+(1-readCol+0.5+lowerCol)*readBuffer[validBand_opt[vband]][lowerCol-startCol];
- if(val_new<=minValue_opt[vband]||val_new>=maxValue_opt[vband]||val_new==srcnodata_opt[vband]){
- readValid=false;
- break;
- }
- }
+ for(int vband=0;vband<bndnodata_opt.size();++vband){
+ val_new=(readCol-0.5-lowerCol)*readBuffer[bndnodata_opt[vband]][upperCol-startCol]+(1-readCol+0.5+lowerCol)*readBuffer[bndnodata_opt[vband]][lowerCol-startCol];
+ if(minValue_opt.size()>vband){
+ if(val_new<=minValue_opt[vband]){
+ readValid=false;
+ break;
+ }
+ }
+ if(maxValue_opt.size()>vband){
+ if(val_new>=maxValue_opt[vband]){
+ readValid=false;
+ break;
+ }
+ }
+ if(srcnodata_opt.size()>vband){
+ if(val_new==srcnodata_opt[vband]){
+ readValid=false;
+ break;
+ }
+ }
+ }
break;
default:
readCol=static_cast<int>(readCol);
- for(int vband=0;vband<validBand_opt.size();++vband){
- val_new=readBuffer[validBand_opt[vband]][readCol-startCol];
- if(val_new<=minValue_opt[vband]||val_new>=maxValue_opt[vband]||val_new==srcnodata_opt[vband]){
- readValid=false;
- break;
- }
- }
+ for(int vband=0;vband<bndnodata_opt.size();++vband){
+ val_new=readBuffer[bndnodata_opt[vband]][readCol-startCol];
+ if(minValue_opt.size()>vband){
+ if(val_new<=minValue_opt[vband]){
+ readValid=false;
+ break;
+ }
+ }
+ if(maxValue_opt.size()>vband){
+ if(val_new>=maxValue_opt[vband]){
+ readValid=false;
+ break;
+ }
+ }
+ if(srcnodata_opt.size()>vband){
+ if(val_new==srcnodata_opt[vband]){
+ readValid=false;
+ break;
+ }
+ }
+ }
break;
- }
+ }
if(readValid){
if(writeValid[ib]){
int iband=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