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 <kempe...@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
_______________________________________________
Pkg-grass-devel mailing list
Pkg-grass-devel@lists.alioth.debian.org
http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-grass-devel
