[pktools] 01/05: Imported Upstream version 2.5.2+20140505
Francesco Lovergine
frankie at moszumanska.debian.org
Wed May 7 07:50:55 UTC 2014
This is an automated email from the git hooks/post-receive script.
frankie pushed a commit to branch master
in repository pktools.
commit 348adb57c8e0cee1125e4aec2f2e675fb1096112
Author: Francesco Paolo Lovergine <frankie at debian.org>
Date: Wed May 7 08:33:24 2014 +0200
Imported Upstream version 2.5.2+20140505
---
ChangeLog | 9 +-
Makefile.am | 24 +
Makefile.in | 108 +-
config.h | 81 --
configure | 29 +-
configure.ac | 27 +-
pktools.pc.in | 11 +
src/algorithms/Filter2d.h | 368 +++++++
src/algorithms/StatFactory.h | 6 +-
src/apps/Makefile.am | 4 +-
src/apps/Makefile.in | 24 +-
src/apps/pkann.cc | 54 +-
src/apps/pkascii2img.cc | 2 +-
src/apps/pkascii2ogr.cc | 2 +-
src/apps/{pkcomposit.cc => pkcomposite.cc} | 42 +-
src/apps/pkcreatect.cc | 2 +-
src/apps/pkcrop.cc | 2 +-
src/apps/pkdiff.cc | 4 +-
src/apps/pkdsm2shadow.cc | 2 +-
src/apps/pkdumpimg.cc | 2 +-
src/apps/pkdumpogr.cc | 2 +-
src/apps/pkdumpogr.h | 2 +-
src/apps/pkeditogr.cc | 2 +-
src/apps/pkegcs.cc | 2 +-
src/apps/pkenhance.cc | 2 +-
src/apps/pkextract.cc | 1481 ++++++++++++----------------
src/apps/pkfillnodata.cc | 2 +-
src/apps/pkfilter.cc | 2 +-
src/apps/pkfilterascii.cc | 2 +-
src/apps/pkfilterdem.cc | 132 ++-
src/apps/pkfsann.cc | 26 +-
src/apps/pkfssvm.cc | 24 +-
src/apps/pkgetmask.cc | 2 +-
src/apps/pkinfo.cc | 50 +-
src/apps/pklas2img.cc | 48 +-
src/apps/pkndvi.cc | 2 +-
src/apps/pkoptsvm.cc | 24 +-
src/apps/pkpolygonize.cc | 2 +-
src/apps/pkreclass.cc | 2 +-
src/apps/pkregann.cc | 4 +-
src/apps/pksetmask.cc | 2 +-
src/apps/pksieve.cc | 2 +-
src/apps/pkstatascii.cc | 4 +-
src/apps/pksvm.cc | 103 +-
src/base/Makefile.am | 2 +-
src/base/Makefile.in | 2 +-
src/base/Vector2d.h | 21 +
src/imageclasses/ImgReaderGdal.cc | 10 +-
src/imageclasses/ImgReaderGdal.h | 2 +-
src/imageclasses/ImgReaderOgr.cc | 32 +-
src/imageclasses/ImgReaderOgr.h | 24 +-
src/imageclasses/ImgWriterOgr.cc | 19 +-
src/imageclasses/ImgWriterOgr.h | 2 +-
src/lasclasses/FileReaderLas.cc | 6 +-
54 files changed, 1631 insertions(+), 1216 deletions(-)
diff --git a/ChangeLog b/ChangeLog
index d0ef20a..16d6767 100755
--- a/ChangeLog
+++ b/ChangeLog
@@ -264,8 +264,11 @@ version 2.5.1
- pksieve
retain nodata in pksieve when mask is set
version 2.5.2
+ - configure script: GDAL>=1.10.0 is required for pkdiff
- programs ported to windows and GUI with Qt
removed underscore for QProcess in Windows
+ - pkinfo
+ distinct long options
- pkclassify_svm -> pksvm
removed underscore for QProcess in Windows
- pkclassify_nn -> pkann
@@ -278,10 +281,12 @@ version 2.5.2
removed underscore for QProcess in Windows
- pkopt_svm -> pkoptsvm
removed underscore for QProcess in Windows
- - pkmosaic -> pkcomposit
+ - pkmosaic -> pkcomposite
name was confusing as also compositing is supported (unlike gdal_merge.py and gdalwarp)
+ resample option similar to pkcrop
- version control for libraries
thanks to suggestion of Francesco Paolo Lovergine
- subdirectory pktools for include headers
thanks to suggestion of Francesco Paolo Lovergine
- -
+ - pklas2img
+ support for compressed point cloud (LAZ) files
diff --git a/Makefile.am b/Makefile.am
index afa0ec6..81d7d64 100755
--- a/Makefile.am
+++ b/Makefile.am
@@ -8,3 +8,27 @@ SUBDIRS = src/base \
src/fileclasses \
$(LASCLASSES_OPT) \
src/apps
+
+## The generated configuration header is installed in its own subdirectory of
+## $(libdir). The reason for this is that the configuration information put
+## into this header file describes the target platform the installed library
+## has been built for. Thus the file must not be installed into a location
+## intended for architecture-independent files, as defined by the Filesystem
+## Hierarchy Standard (FHS).
+## The nodist_ prefix instructs Automake to not generate rules for including
+## the listed files in the distribution on 'make dist'. Files that are listed
+## in _HEADERS variables are normally included in the distribution, but the
+## configuration header file is generated at configure time and should not be
+## shipped with the source tarball.
+pktools_libincludedir = $(libdir)/pktools/include
+nodist_pktools_libinclude_HEADERS = config.h
+
+
+## Install the generated pkg-config file (.pc) into the expected location for
+## architecture-dependent package configuration information. Occasionally,
+## pkg-config files are also used for architecture-independent data packages,
+## in which case the correct install location would be $(datadir)/pkgconfig.
+pkgconfigdir = $(libdir)/pkgconfig
+pkgconfig_DATA = pktools.pc
+
+#bin_SCRIPTS = pktools-config
diff --git a/Makefile.in b/Makefile.in
index 5bd6c1b..a2b41a4 100644
--- a/Makefile.in
+++ b/Makefile.in
@@ -14,6 +14,8 @@
# PARTICULAR PURPOSE.
@SET_MAKE@
+
+
VPATH = @srcdir@
am__make_dryrun = \
{ \
@@ -53,8 +55,9 @@ host_triplet = @host@
subdir = .
DIST_COMMON = README $(am__configure_deps) $(srcdir)/Makefile.am \
$(srcdir)/Makefile.in $(srcdir)/config.h.in \
- $(top_srcdir)/configure AUTHORS COPYING ChangeLog INSTALL NEWS \
- config.guess config.sub depcomp install-sh ltmain.sh missing
+ $(srcdir)/pktools.pc.in $(top_srcdir)/configure AUTHORS \
+ COPYING ChangeLog INSTALL NEWS config.guess config.sub depcomp \
+ install-sh ltmain.sh missing
ACLOCAL_M4 = $(top_srcdir)/aclocal.m4
am__aclocal_m4_deps = $(top_srcdir)/m4/ax_lib_gdal.m4 \
$(top_srcdir)/m4/libtool.m4 $(top_srcdir)/m4/ltoptions.m4 \
@@ -66,7 +69,7 @@ am__CONFIG_DISTCLEAN_FILES = config.status config.cache config.log \
configure.lineno config.status.lineno
mkinstalldirs = $(install_sh) -d
CONFIG_HEADER = config.h
-CONFIG_CLEAN_FILES =
+CONFIG_CLEAN_FILES = pktools.pc
CONFIG_CLEAN_VPATH_FILES =
SOURCES =
DIST_SOURCES =
@@ -82,6 +85,37 @@ am__can_run_installinfo = \
n|no|NO) false;; \
*) (install-info --version) >/dev/null 2>&1;; \
esac
+am__vpath_adj_setup = srcdirstrip=`echo "$(srcdir)" | sed 's|.|.|g'`;
+am__vpath_adj = case $$p in \
+ $(srcdir)/*) f=`echo "$$p" | sed "s|^$$srcdirstrip/||"`;; \
+ *) f=$$p;; \
+ esac;
+am__strip_dir = f=`echo $$p | sed -e 's|^.*/||'`;
+am__install_max = 40
+am__nobase_strip_setup = \
+ srcdirstrip=`echo "$(srcdir)" | sed 's/[].[^$$\\*|]/\\\\&/g'`
+am__nobase_strip = \
+ for p in $$list; do echo "$$p"; done | sed -e "s|$$srcdirstrip/||"
+am__nobase_list = $(am__nobase_strip_setup); \
+ for p in $$list; do echo "$$p $$p"; done | \
+ sed "s| $$srcdirstrip/| |;"' / .*\//!s/ .*/ ./; s,\( .*\)/[^/]*$$,\1,' | \
+ $(AWK) 'BEGIN { files["."] = "" } { files[$$2] = files[$$2] " " $$1; \
+ if (++n[$$2] == $(am__install_max)) \
+ { print $$2, files[$$2]; n[$$2] = 0; files[$$2] = "" } } \
+ END { for (dir in files) print dir, files[dir] }'
+am__base_list = \
+ sed '$$!N;$$!N;$$!N;$$!N;$$!N;$$!N;$$!N;s/\n/ /g' | \
+ sed '$$!N;$$!N;$$!N;$$!N;s/\n/ /g'
+am__uninstall_files_from_dir = { \
+ test -z "$$files" \
+ || { test ! -d "$$dir" && test ! -f "$$dir" && test ! -r "$$dir"; } \
+ || { echo " ( cd '$$dir' && rm -f" $$files ")"; \
+ $(am__cd) "$$dir" && rm -f $$files; }; \
+ }
+am__installdirs = "$(DESTDIR)$(pkgconfigdir)" \
+ "$(DESTDIR)$(pktools_libincludedir)"
+DATA = $(pkgconfig_DATA)
+HEADERS = $(nodist_pktools_libinclude_HEADERS)
RECURSIVE_CLEAN_TARGETS = mostlyclean-recursive clean-recursive \
distclean-recursive maintainer-clean-recursive
AM_RECURSIVE_TARGETS = $(RECURSIVE_TARGETS:-recursive=) \
@@ -282,6 +316,10 @@ SUBDIRS = src/base \
$(LASCLASSES_OPT) \
src/apps
+pktools_libincludedir = $(libdir)/pktools/include
+nodist_pktools_libinclude_HEADERS = config.h
+pkgconfigdir = $(libdir)/pkgconfig
+pkgconfig_DATA = pktools.pc
all: config.h
$(MAKE) $(AM_MAKEFLAGS) all-recursive
@@ -335,6 +373,8 @@ $(srcdir)/config.h.in: $(am__configure_deps)
distclean-hdr:
-rm -f config.h stamp-h1
+pktools.pc: $(top_builddir)/config.status $(srcdir)/pktools.pc.in
+ cd $(top_builddir) && $(SHELL) ./config.status $@
mostlyclean-libtool:
-rm -f *.lo
@@ -344,6 +384,48 @@ clean-libtool:
distclean-libtool:
-rm -f libtool config.lt
+install-pkgconfigDATA: $(pkgconfig_DATA)
+ @$(NORMAL_INSTALL)
+ @list='$(pkgconfig_DATA)'; test -n "$(pkgconfigdir)" || list=; \
+ if test -n "$$list"; then \
+ echo " $(MKDIR_P) '$(DESTDIR)$(pkgconfigdir)'"; \
+ $(MKDIR_P) "$(DESTDIR)$(pkgconfigdir)" || exit 1; \
+ fi; \
+ for p in $$list; do \
+ if test -f "$$p"; then d=; else d="$(srcdir)/"; fi; \
+ echo "$$d$$p"; \
+ done | $(am__base_list) | \
+ while read files; do \
+ echo " $(INSTALL_DATA) $$files '$(DESTDIR)$(pkgconfigdir)'"; \
+ $(INSTALL_DATA) $$files "$(DESTDIR)$(pkgconfigdir)" || exit $$?; \
+ done
+
+uninstall-pkgconfigDATA:
+ @$(NORMAL_UNINSTALL)
+ @list='$(pkgconfig_DATA)'; test -n "$(pkgconfigdir)" || list=; \
+ files=`for p in $$list; do echo $$p; done | sed -e 's|^.*/||'`; \
+ dir='$(DESTDIR)$(pkgconfigdir)'; $(am__uninstall_files_from_dir)
+install-nodist_pktools_libincludeHEADERS: $(nodist_pktools_libinclude_HEADERS)
+ @$(NORMAL_INSTALL)
+ @list='$(nodist_pktools_libinclude_HEADERS)'; test -n "$(pktools_libincludedir)" || list=; \
+ if test -n "$$list"; then \
+ echo " $(MKDIR_P) '$(DESTDIR)$(pktools_libincludedir)'"; \
+ $(MKDIR_P) "$(DESTDIR)$(pktools_libincludedir)" || exit 1; \
+ fi; \
+ for p in $$list; do \
+ if test -f "$$p"; then d=; else d="$(srcdir)/"; fi; \
+ echo "$$d$$p"; \
+ done | $(am__base_list) | \
+ while read files; do \
+ echo " $(INSTALL_HEADER) $$files '$(DESTDIR)$(pktools_libincludedir)'"; \
+ $(INSTALL_HEADER) $$files "$(DESTDIR)$(pktools_libincludedir)" || exit $$?; \
+ done
+
+uninstall-nodist_pktools_libincludeHEADERS:
+ @$(NORMAL_UNINSTALL)
+ @list='$(nodist_pktools_libinclude_HEADERS)'; test -n "$(pktools_libincludedir)" || list=; \
+ files=`for p in $$list; do echo $$p; done | sed -e 's|^.*/||'`; \
+ dir='$(DESTDIR)$(pktools_libincludedir)'; $(am__uninstall_files_from_dir)
# This directory's subdirectories are mostly independent; you can cd
# into them and run `make' without going through this Makefile.
@@ -669,9 +751,12 @@ distcleancheck: distclean
exit 1; } >&2
check-am: all-am
check: check-recursive
-all-am: Makefile config.h
+all-am: Makefile $(DATA) $(HEADERS) config.h
installdirs: installdirs-recursive
installdirs-am:
+ for dir in "$(DESTDIR)$(pkgconfigdir)" "$(DESTDIR)$(pktools_libincludedir)"; do \
+ test -z "$$dir" || $(MKDIR_P) "$$dir"; \
+ done
install: install-recursive
install-exec: install-exec-recursive
install-data: install-data-recursive
@@ -724,7 +809,8 @@ info: info-recursive
info-am:
-install-data-am:
+install-data-am: install-nodist_pktools_libincludeHEADERS \
+ install-pkgconfigDATA
install-dvi: install-dvi-recursive
@@ -770,7 +856,8 @@ ps: ps-recursive
ps-am:
-uninstall-am:
+uninstall-am: uninstall-nodist_pktools_libincludeHEADERS \
+ uninstall-pkgconfigDATA
.MAKE: $(RECURSIVE_CLEAN_TARGETS) $(RECURSIVE_TARGETS) all \
ctags-recursive install-am install-strip tags-recursive
@@ -785,12 +872,17 @@ uninstall-am:
install install-am install-data install-data-am install-dvi \
install-dvi-am install-exec install-exec-am install-html \
install-html-am install-info install-info-am install-man \
- install-pdf install-pdf-am install-ps install-ps-am \
+ install-nodist_pktools_libincludeHEADERS install-pdf \
+ install-pdf-am install-pkgconfigDATA install-ps install-ps-am \
install-strip installcheck installcheck-am installdirs \
installdirs-am maintainer-clean maintainer-clean-generic \
mostlyclean mostlyclean-generic mostlyclean-libtool pdf pdf-am \
- ps ps-am tags tags-recursive uninstall uninstall-am
+ ps ps-am tags tags-recursive uninstall uninstall-am \
+ uninstall-nodist_pktools_libincludeHEADERS \
+ uninstall-pkgconfigDATA
+
+#bin_SCRIPTS = pktools-config
# Tell versions [3.59,3.63) of GNU make to not export all variables.
# Otherwise a system limit (for SysV at least) may be exceeded.
diff --git a/config.h b/config.h
deleted file mode 100644
index 975dcb0..0000000
--- a/config.h
+++ /dev/null
@@ -1,81 +0,0 @@
-/* config.h. Generated from config.h.in by configure. */
-/* config.h.in. Generated from configure.ac by autoheader. */
-
-/* Define to 1 if you have the <dlfcn.h> header file. */
-#define HAVE_DLFCN_H 1
-
-/* Define to 1 if GDAL library are available */
-#define HAVE_GDAL 1
-
-/* Define to 1 if you have the <gdal.h> header file. */
-/* #undef HAVE_GDAL_H */
-
-/* Define to 1 if GDAL library includes OGR support */
-#define HAVE_GDAL_OGR 1
-
-/* Define to 1 if you have the <inttypes.h> header file. */
-#define HAVE_INTTYPES_H 1
-
-/* Define to 1 if you have the <iostream> header file. */
-#define HAVE_IOSTREAM 1
-
-/* Define to 1 if you have the <memory.h> header file. */
-#define HAVE_MEMORY_H 1
-
-/* Define to 1 if you have the <stdint.h> header file. */
-#define HAVE_STDINT_H 1
-
-/* Define to 1 if you have the <stdlib.h> header file. */
-#define HAVE_STDLIB_H 1
-
-/* Define to 1 if you have the <string> header file. */
-#define HAVE_STRING 1
-
-/* Define to 1 if you have the <strings.h> header file. */
-#define HAVE_STRINGS_H 1
-
-/* Define to 1 if you have the <string.h> header file. */
-#define HAVE_STRING_H 1
-
-/* Define to 1 if you have the <sys/stat.h> header file. */
-#define HAVE_SYS_STAT_H 1
-
-/* Define to 1 if you have the <sys/types.h> header file. */
-#define HAVE_SYS_TYPES_H 1
-
-/* Define to 1 if you have the <unistd.h> header file. */
-#define HAVE_UNISTD_H 1
-
-/* Define to the sub-directory in which libtool stores uninstalled libraries.
- */
-#define LT_OBJDIR ".libs/"
-
-/* Name of package */
-#define PACKAGE "pktools"
-
-/* Define to the address where bug reports for this package should be sent. */
-#define PACKAGE_BUGREPORT "kempenep at gmail.com"
-
-/* Define to the full name of this package. */
-#define PACKAGE_NAME "pktools"
-
-/* Define to the full name and version of this package. */
-#define PACKAGE_STRING "pktools 2.5.2"
-
-/* Define to the one symbol short name of this package. */
-#define PACKAGE_TARNAME "pktools"
-
-/* Define to the home page for this package. */
-#define PACKAGE_URL ""
-
-/* Define to the version of this package. */
-#define PACKAGE_VERSION "2.5.2"
-
-/* Define to 1 if you have the ANSI C header files. */
-#define STDC_HEADERS 1
-
-/* Version number of package */
-#define VERSION "2.5.2"
-
-/* Define to `unsigned int' if <sys/types.h> does not define. */
-/* #undef size_t */
diff --git a/configure b/configure
index 05b6c09..e8cb6b9 100755
--- a/configure
+++ b/configure
@@ -3366,16 +3366,6 @@ GDAL_FOUND="no"
if test ! -z "$GDAL_CFLAGS" -a ! -z "$GDAL_LDFLAGS"; then
GDAL_FOUND="yes"
fi
-{ $as_echo "$as_me:${as_lineno-$LINENO}: WARNING: $GDAL_VERSION" >&5
-$as_echo "$as_me: WARNING: $GDAL_VERSION" >&2;}
-{ $as_echo "$as_me:${as_lineno-$LINENO}: WARNING: $GDAL_CFLAGS" >&5
-$as_echo "$as_me: WARNING: $GDAL_CFLAGS" >&2;}
-{ $as_echo "$as_me:${as_lineno-$LINENO}: WARNING: $GDAL_LDFLAGS" >&5
-$as_echo "$as_me: WARNING: $GDAL_LDFLAGS" >&2;}
-{ $as_echo "$as_me:${as_lineno-$LINENO}: WARNING: $GDAL_DEP_LDFLAGS" >&5
-$as_echo "$as_me: WARNING: $GDAL_DEP_LDFLAGS" >&2;}
-{ $as_echo "$as_me:${as_lineno-$LINENO}: WARNING: $GDAL_OGR_ENABLED" >&5
-$as_echo "$as_me: WARNING: $GDAL_OGR_ENABLED" >&2;}
# check for C++ compiler and the library compiler
ac_ext=cpp
@@ -18764,7 +18754,7 @@ $as_echo "$GDAL_OGR_ENABLED" >&6; }
fi
- gdal_version_req=1.8.0
+ gdal_version_req=1.10.0
if test "$found_gdal" = "yes" -a -n "$gdal_version_req"; then
{ $as_echo "$as_me:${as_lineno-$LINENO}: checking if GDAL version is >= $gdal_version_req" >&5
@@ -19029,8 +19019,6 @@ else
USE_FANN_FALSE=
fi
-{ $as_echo "$as_me:${as_lineno-$LINENO}: WARNING: \"USE_FANN=true\"" >&5
-$as_echo "$as_me: WARNING: \"USE_FANN=true\"" >&2;}
else
@@ -19042,8 +19030,6 @@ else
USE_FANN_FALSE=
fi
-{ $as_echo "$as_me:${as_lineno-$LINENO}: WARNING: \"USE_FANN=false\"" >&5
-$as_echo "$as_me: WARNING: \"USE_FANN=false\"" >&2;}
fi
@@ -19065,8 +19051,6 @@ else
USE_LAS_FALSE=
fi
-{ $as_echo "$as_me:${as_lineno-$LINENO}: WARNING: \"USE_LAS=true\"" >&5
-$as_echo "$as_me: WARNING: \"USE_LAS=true\"" >&2;}
LASCLASSES_OPT="src/lasclasses"
@@ -19080,8 +19064,6 @@ else
USE_LAS_FALSE=
fi
-{ $as_echo "$as_me:${as_lineno-$LINENO}: WARNING: \"USE_LAS=false\"" >&5
-$as_echo "$as_me: WARNING: \"USE_LAS=false\"" >&2;}
fi
@@ -19172,8 +19154,6 @@ else
USE_NLOPT_FALSE=
fi
-{ $as_echo "$as_me:${as_lineno-$LINENO}: WARNING: \"USE_NLOPT=true\"" >&5
-$as_echo "$as_me: WARNING: \"USE_NLOPT=true\"" >&2;}
else
@@ -19185,8 +19165,6 @@ else
USE_NLOPT_FALSE=
fi
-{ $as_echo "$as_me:${as_lineno-$LINENO}: WARNING: \"USE_NLOPT=false\"" >&5
-$as_echo "$as_me: WARNING: \"USE_NLOPT=false\"" >&2;}
fi
@@ -19447,12 +19425,14 @@ PKTOOLS_SO_VERSION=1:0:0
# files to generate via autotools (.am or .in source files)
ac_config_headers="$ac_config_headers config.h"
+#frankie says: "You should not distribute the config.h file, but the config.h.in file only".
+#AC_CONFIG_HEADERS([config.h])
if test -z "$USE_LAS_TRUE"; then :
ac_config_files="$ac_config_files src/lasclasses/Makefile"
fi
-ac_config_files="$ac_config_files Makefile src/base/Makefile src/algorithms/Makefile src/imageclasses/Makefile src/fileclasses/Makefile src/apps/Makefile"
+ac_config_files="$ac_config_files Makefile src/base/Makefile src/algorithms/Makefile src/imageclasses/Makefile src/fileclasses/Makefile src/apps/Makefile pktools.pc:pktools.pc.in"
# generate the final Makefile etc.
cat >confcache <<\_ACEOF
@@ -20675,6 +20655,7 @@ do
"src/imageclasses/Makefile") CONFIG_FILES="$CONFIG_FILES src/imageclasses/Makefile" ;;
"src/fileclasses/Makefile") CONFIG_FILES="$CONFIG_FILES src/fileclasses/Makefile" ;;
"src/apps/Makefile") CONFIG_FILES="$CONFIG_FILES src/apps/Makefile" ;;
+ "pktools.pc") CONFIG_FILES="$CONFIG_FILES pktools.pc:pktools.pc.in" ;;
*) as_fn_error $? "invalid argument: \`$ac_config_target'" "$LINENO" 5;;
esac
diff --git a/configure.ac b/configure.ac
index fae3f47..30ec14a 100644
--- a/configure.ac
+++ b/configure.ac
@@ -8,11 +8,11 @@ GDAL_FOUND="no"
if test ! -z "$GDAL_CFLAGS" -a ! -z "$GDAL_LDFLAGS"; then
GDAL_FOUND="yes"
fi
-AC_MSG_WARN([$GDAL_VERSION])
-AC_MSG_WARN([$GDAL_CFLAGS])
-AC_MSG_WARN([$GDAL_LDFLAGS])
-AC_MSG_WARN([$GDAL_DEP_LDFLAGS])
-AC_MSG_WARN([$GDAL_OGR_ENABLED])
+dnl AC_MSG_WARN([$GDAL_VERSION])
+dnl AC_MSG_WARN([$GDAL_CFLAGS])
+dnl AC_MSG_WARN([$GDAL_LDFLAGS])
+dnl AC_MSG_WARN([$GDAL_DEP_LDFLAGS])
+dnl AC_MSG_WARN([$GDAL_OGR_ENABLED])
# check for C++ compiler and the library compiler
AC_PROG_CXX
@@ -21,7 +21,7 @@ LT_INIT
# check if the source folder is correct
AC_CONFIG_SRCDIR([src/apps/pkinfo.cc])
-AX_LIB_GDAL([1.8.0]) dnl uncomment if gdal version 1.8.0 is required
+AX_LIB_GDAL([1.10.0]) dnl uncomment if gdal version 1.10.0 is required
AC_CHECK_HEADERS([gdal.h])
@@ -32,10 +32,10 @@ AS_HELP_STRING([--enable-fann], [Enable feature fann]))
AS_IF([test "x$enable_fann" = "xyes"], [
PKG_CHECK_MODULES(FANN, fann >= 2.1.0, [HAVE_FANN=1], [HAVE_FANN=0])
AM_CONDITIONAL([USE_FANN], [test "$HAVE_FANN" -eq 1])
-AC_MSG_WARN("USE_FANN=true")
+dnl AC_MSG_WARN("USE_FANN=true")
], [
AM_CONDITIONAL(USE_FANN, false)
-AC_MSG_WARN("USE_FANN=false")
+dnl AC_MSG_WARN("USE_FANN=false")
])
dnl AC_ARG_ENABLE([svm],
@@ -51,12 +51,12 @@ AS_HELP_STRING([--enable-las], [Enable feature las]))
AS_IF([test "x$enable_las" = "xyes"], [
AM_CONDITIONAL(USE_LAS, true)
-AC_MSG_WARN("USE_LAS=true")
+dnl AC_MSG_WARN("USE_LAS=true")
LASCLASSES_OPT="src/lasclasses"
AC_SUBST([LASCLASSES_OPT])
], [
AM_CONDITIONAL(USE_LAS, false)
-AC_MSG_WARN("USE_LAS=false")
+dnl AC_MSG_WARN("USE_LAS=false")
])
AC_ARG_ENABLE([nlopt],
@@ -65,10 +65,10 @@ AS_HELP_STRING([--enable-nlopt], [Enable feature nlopt]))
AS_IF([test "x$enable_nlopt" = "xyes"], [
PKG_CHECK_MODULES(NLOPT, nlopt >= 2.1.0, [HAVE_NLOPT=1], [HAVE_NLOPT=0])
AM_CONDITIONAL([USE_NLOPT], [test "$HAVE_NLOPT" -eq 1])
-AC_MSG_WARN("USE_NLOPT=true")
+dnl AC_MSG_WARN("USE_NLOPT=true")
], [
AM_CONDITIONAL(USE_NLOPT, false)
-AC_MSG_WARN("USE_NLOPT=false")
+dnl AC_MSG_WARN("USE_NLOPT=false")
])
@@ -100,6 +100,8 @@ AC_SUBST([PKTOOLS_SO_VERSION], [1:0:0])
# files to generate via autotools (.am or .in source files)
AC_CONFIG_HEADERS([config.h])
+#frankie says: "You should not distribute the config.h file, but the config.h.in file only".
+#AC_CONFIG_HEADERS([config.h])
AM_COND_IF([USE_LAS],
[AC_CONFIG_FILES([
@@ -111,6 +113,7 @@ AC_CONFIG_FILES([
src/imageclasses/Makefile
src/fileclasses/Makefile
src/apps/Makefile
+ pktools.pc:pktools.pc.in
])
# generate the final Makefile etc.
AC_OUTPUT
diff --git a/pktools.pc.in b/pktools.pc.in
new file mode 100644
index 0000000..0731529
--- /dev/null
+++ b/pktools.pc.in
@@ -0,0 +1,11 @@
+prefix=@prefix@
+exec_prefix=@exec_prefix@
+libdir=@libdir@
+includedir=@includedir@
+
+Name: pktools
+Description: API library for pktools
+Requires: gdal gsl
+Version: @PACKAGE_VERSION@
+Libs: -L${libdir} -lpktools
+Cflags: -I${includedir}/pktools
\ No newline at end of file
diff --git a/src/algorithms/Filter2d.h b/src/algorithms/Filter2d.h
index 8d3717c..54aeee0 100644
--- a/src/algorithms/Filter2d.h
+++ b/src/algorithms/Filter2d.h
@@ -116,6 +116,10 @@ public:
void var(const std::string& inputFilename, const std::string& outputFilename, int dim, bool disc=false);
void morphology(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, int dimX, int dimY, const std::vector<double> &angle, bool disc=false);
template<class T> unsigned long int morphology(const Vector2d<T>& input, Vector2d<T>& output, const std::string& method, int dimX, int dimY, bool disc=false, double hThreshold=0);
+ template<class T> unsigned long int dsm2dtm_nwse(const Vector2d<T>& inputDSM, Vector2d<T>& outputMask, double hThreshold, int nlimit, int dim=3);
+ template<class T> unsigned long int dsm2dtm_nesw(const Vector2d<T>& inputDSM, Vector2d<T>& outputMask, double hThreshold, int nlimit, int dim=3);
+ template<class T> unsigned long int dsm2dtm_senw(const Vector2d<T>& inputDSM, Vector2d<T>& outputMask, double hThreshold, int nlimit, int dim=3);
+ template<class T> unsigned long int dsm2dtm_swne(const Vector2d<T>& inputDSM, Vector2d<T>& outputMask, double hThreshold, int nlimit, int dim=3);
template<class T> void shadowDsm(const Vector2d<T>& input, Vector2d<T>& output, double sza, double saa, double pixelSize, short shadowFlag=1);
void shadowDsm(const ImgReaderGdal& input, ImgWriterGdal& output, double sza, double saa, double pixelSize, short shadowFlag=1);
void dwt_texture(const std::string& inputFilename, const std::string& outputFilename,int dim, int scale, int down=1, int iband=0, bool verbose=false);
@@ -747,6 +751,370 @@ template<class T> unsigned long int Filter2d::morphology(const Vector2d<T>& inpu
return nchange;
}
+ template<class T> unsigned long int Filter2d::dsm2dtm_nwse(const Vector2d<T>& inputDSM, Vector2d<T>& outputMask, double hThreshold, int nlimit, int dim)
+{
+ const char* pszMessage;
+ void* pProgressArg=NULL;
+ GDALProgressFunc pfnProgress=GDALTermProgress;
+ double progress=0;
+ pfnProgress(progress,pszMessage,pProgressArg);
+
+ Vector2d<T> tmpDSM(inputDSM);
+ double noDataValue=0;
+ if(m_noDataValues.size())
+ noDataValue=m_noDataValues[0];
+
+ unsigned long int nchange=0;
+ int dimX=dim;
+ int dimY=dim;
+ assert(dimX);
+ assert(dimY);
+ statfactory::StatFactory stat;
+ Vector2d<T> inBuffer(dimY,inputDSM.nCols());
+ if(outputMask.size()!=inputDSM.nRows())
+ outputMask.resize(inputDSM.nRows());
+ int indexI=0;
+ int indexJ=0;
+ //initialize last half of inBuffer
+ for(int j=-(dimY-1)/2;j<=dimY/2;++j){
+ for(int i=0;i<inputDSM.nCols();++i)
+ inBuffer[indexJ][i]=tmpDSM[abs(j)][i];
+ ++indexJ;
+ }
+ for(int y=0;y<tmpDSM.nRows();++y){
+ if(y){//inBuffer already initialized for y=0
+ //erase first line from inBuffer
+ inBuffer.erase(inBuffer.begin());
+ //read extra line and push back to inBuffer if not out of bounds
+ if(y+dimY/2<tmpDSM.nRows()){
+ //allocate buffer
+ inBuffer.push_back(inBuffer.back());
+ for(int i=0;i<tmpDSM.nCols();++i)
+ inBuffer[inBuffer.size()-1][i]=tmpDSM[y+dimY/2][i];
+ }
+ else{
+ int over=y+dimY/2-tmpDSM.nRows();
+ int index=(inBuffer.size()-1)-over;
+ assert(index>=0);
+ assert(index<inBuffer.size());
+ inBuffer.push_back(inBuffer[index]);
+ }
+ }
+ for(int x=0;x<tmpDSM.nCols();++x){
+ double centerValue=inBuffer[(dimY-1)/2][x];
+ short nmasked=0;
+ std::vector<T> neighbors;
+ for(int j=-(dimY-1)/2;j<=dimY/2;++j){
+ for(int i=-(dimX-1)/2;i<=dimX/2;++i){
+ indexI=x+i;
+ //check if out of bounds
+ if(indexI<0)
+ indexI=-indexI;
+ else if(indexI>=tmpDSM.nCols())
+ indexI=tmpDSM.nCols()-i;
+ if(y+j<0)
+ indexJ=-j;
+ else if(y+j>=tmpDSM.nRows())
+ indexJ=(dimY>2) ? (dimY-1)/2-j : 0;
+ else
+ indexJ=(dimY-1)/2+j;
+ double difference=(centerValue-inBuffer[indexJ][indexI]);
+ if(i||j)//skip centerValue
+ neighbors.push_back(inBuffer[indexJ][indexI]);
+ if(difference>hThreshold)
+ ++nmasked;
+ }
+ }
+ if(nmasked<nlimit){
+ ++nchange;
+ //reset pixel in outputMask
+ outputMask[y][x]=0;
+ }
+ else{
+ //reset pixel height in tmpDSM
+ inBuffer[(dimY-1)/2][x]=stat.mymin(neighbors);
+ }
+ }
+ progress=(1.0+y);
+ progress/=outputMask.nRows();
+ pfnProgress(progress,pszMessage,pProgressArg);
+ }
+ return nchange;
+}
+
+ template<class T> unsigned long int Filter2d::dsm2dtm_nesw(const Vector2d<T>& inputDSM, Vector2d<T>& outputMask, double hThreshold, int nlimit, int dim)
+{
+ const char* pszMessage;
+ void* pProgressArg=NULL;
+ GDALProgressFunc pfnProgress=GDALTermProgress;
+ double progress=0;
+ pfnProgress(progress,pszMessage,pProgressArg);
+
+ Vector2d<T> tmpDSM(inputDSM);
+ double noDataValue=0;
+ if(m_noDataValues.size())
+ noDataValue=m_noDataValues[0];
+
+ unsigned long int nchange=0;
+ int dimX=dim;
+ int dimY=dim;
+ assert(dimX);
+ assert(dimY);
+ statfactory::StatFactory stat;
+ Vector2d<T> inBuffer(dimY,inputDSM.nCols());
+ if(outputMask.size()!=inputDSM.nRows())
+ outputMask.resize(inputDSM.nRows());
+ int indexI=0;
+ int indexJ=0;
+ //initialize last half of inBuffer
+ for(int j=-(dimY-1)/2;j<=dimY/2;++j){
+ for(int i=0;i<inputDSM.nCols();++i)
+ inBuffer[indexJ][i]=tmpDSM[abs(j)][i];
+ ++indexJ;
+ }
+ for(int y=0;y<tmpDSM.nRows();++y){
+ if(y){//inBuffer already initialized for y=0
+ //erase first line from inBuffer
+ inBuffer.erase(inBuffer.begin());
+ //read extra line and push back to inBuffer if not out of bounds
+ if(y+dimY/2<tmpDSM.nRows()){
+ //allocate buffer
+ inBuffer.push_back(inBuffer.back());
+ for(int i=0;i<tmpDSM.nCols();++i)
+ inBuffer[inBuffer.size()-1][i]=tmpDSM[y+dimY/2][i];
+ }
+ else{
+ int over=y+dimY/2-tmpDSM.nRows();
+ int index=(inBuffer.size()-1)-over;
+ assert(index>=0);
+ assert(index<inBuffer.size());
+ inBuffer.push_back(inBuffer[index]);
+ }
+ }
+ for(int x=tmpDSM.nCols()-1;x>=0;--x){
+ double centerValue=inBuffer[(dimY-1)/2][x];
+ short nmasked=0;
+ std::vector<T> neighbors;
+ for(int j=-(dimY-1)/2;j<=dimY/2;++j){
+ for(int i=-(dimX-1)/2;i<=dimX/2;++i){
+ indexI=x+i;
+ //check if out of bounds
+ if(indexI<0)
+ indexI=-indexI;
+ else if(indexI>=tmpDSM.nCols())
+ indexI=tmpDSM.nCols()-i;
+ if(y+j<0)
+ indexJ=-j;
+ else if(y+j>=tmpDSM.nRows())
+ indexJ=(dimY>2) ? (dimY-1)/2-j : 0;
+ else
+ indexJ=(dimY-1)/2+j;
+ double difference=(centerValue-inBuffer[indexJ][indexI]);
+ if(i||j)//skip centerValue
+ neighbors.push_back(inBuffer[indexJ][indexI]);
+ if(difference>hThreshold)
+ ++nmasked;
+ }
+ }
+ if(nmasked<nlimit){
+ ++nchange;
+ //reset pixel in outputMask
+ outputMask[y][x]=0;
+ }
+ else{
+ //reset pixel height in tmpDSM
+ inBuffer[(dimY-1)/2][x]=stat.mymin(neighbors);
+ }
+ }
+ progress=(1.0+y);
+ progress/=outputMask.nRows();
+ pfnProgress(progress,pszMessage,pProgressArg);
+ }
+ return nchange;
+}
+
+ template<class T> unsigned long int Filter2d::dsm2dtm_senw(const Vector2d<T>& inputDSM, Vector2d<T>& outputMask, double hThreshold, int nlimit, int dim)
+{
+ const char* pszMessage;
+ void* pProgressArg=NULL;
+ GDALProgressFunc pfnProgress=GDALTermProgress;
+ double progress=0;
+ pfnProgress(progress,pszMessage,pProgressArg);
+
+ Vector2d<T> tmpDSM(inputDSM);
+ double noDataValue=0;
+ if(m_noDataValues.size())
+ noDataValue=m_noDataValues[0];
+
+ unsigned long int nchange=0;
+ int dimX=dim;
+ int dimY=dim;
+ assert(dimX);
+ assert(dimY);
+ statfactory::StatFactory stat;
+ Vector2d<T> inBuffer(dimY,inputDSM.nCols());
+ if(outputMask.size()!=inputDSM.nRows())
+ outputMask.resize(inputDSM.nRows());
+ int indexI=0;
+ int indexJ=0;
+ //initialize last half of inBuffer
+ for(int j=-(dimY-1)/2;j<=dimY/2;++j){
+ for(int i=0;i<inputDSM.nCols();++i)
+ inBuffer[indexJ][i]=tmpDSM[abs(j)][i];
+ ++indexJ;
+ }
+ for(int y=tmpDSM.nRows()-1;y>=0;--y){
+ if(y){//inBuffer already initialized for y=0
+ //erase first line from inBuffer
+ inBuffer.erase(inBuffer.begin());
+ //read extra line and push back to inBuffer if not out of bounds
+ if(y+dimY/2<tmpDSM.nRows()){
+ //allocate buffer
+ inBuffer.push_back(inBuffer.back());
+ for(int i=0;i<tmpDSM.nCols();++i)
+ inBuffer[inBuffer.size()-1][i]=tmpDSM[y+dimY/2][i];
+ }
+ else{
+ int over=y+dimY/2-tmpDSM.nRows();
+ int index=(inBuffer.size()-1)-over;
+ assert(index>=0);
+ assert(index<inBuffer.size());
+ inBuffer.push_back(inBuffer[index]);
+ }
+ }
+ for(int x=tmpDSM.nCols()-1;x>=0;--x){
+ double centerValue=inBuffer[(dimY-1)/2][x];
+ short nmasked=0;
+ std::vector<T> neighbors;
+ for(int j=-(dimY-1)/2;j<=dimY/2;++j){
+ for(int i=-(dimX-1)/2;i<=dimX/2;++i){
+ indexI=x+i;
+ //check if out of bounds
+ if(indexI<0)
+ indexI=-indexI;
+ else if(indexI>=tmpDSM.nCols())
+ indexI=tmpDSM.nCols()-i;
+ if(y+j<0)
+ indexJ=-j;
+ else if(y+j>=tmpDSM.nRows())
+ indexJ=(dimY>2) ? (dimY-1)/2-j : 0;
+ else
+ indexJ=(dimY-1)/2+j;
+ double difference=(centerValue-inBuffer[indexJ][indexI]);
+ if(i||j)//skip centerValue
+ neighbors.push_back(inBuffer[indexJ][indexI]);
+ if(difference>hThreshold)
+ ++nmasked;
+ }
+ }
+ if(nmasked<nlimit){
+ ++nchange;
+ //reset pixel in outputMask
+ outputMask[y][x]=0;
+ }
+ else{
+ //reset pixel height in tmpDSM
+ inBuffer[(dimY-1)/2][x]=stat.mymin(neighbors);
+ }
+ }
+ progress=(1.0+y);
+ progress/=outputMask.nRows();
+ pfnProgress(progress,pszMessage,pProgressArg);
+ }
+ return nchange;
+}
+
+ template<class T> unsigned long int Filter2d::dsm2dtm_swne(const Vector2d<T>& inputDSM, Vector2d<T>& outputMask, double hThreshold, int nlimit, int dim)
+{
+ const char* pszMessage;
+ void* pProgressArg=NULL;
+ GDALProgressFunc pfnProgress=GDALTermProgress;
+ double progress=0;
+ pfnProgress(progress,pszMessage,pProgressArg);
+
+ Vector2d<T> tmpDSM(inputDSM);
+ double noDataValue=0;
+ if(m_noDataValues.size())
+ noDataValue=m_noDataValues[0];
+
+ unsigned long int nchange=0;
+ int dimX=dim;
+ int dimY=dim;
+ assert(dimX);
+ assert(dimY);
+ statfactory::StatFactory stat;
+ Vector2d<T> inBuffer(dimY,inputDSM.nCols());
+ if(outputMask.size()!=inputDSM.nRows())
+ outputMask.resize(inputDSM.nRows());
+ int indexI=0;
+ int indexJ=0;
+ //initialize last half of inBuffer
+ for(int j=-(dimY-1)/2;j<=dimY/2;++j){
+ for(int i=0;i<inputDSM.nCols();++i)
+ inBuffer[indexJ][i]=tmpDSM[abs(j)][i];
+ ++indexJ;
+ }
+ for(int y=tmpDSM.nRows()-1;y>=0;--y){
+ if(y){//inBuffer already initialized for y=0
+ //erase first line from inBuffer
+ inBuffer.erase(inBuffer.begin());
+ //read extra line and push back to inBuffer if not out of bounds
+ if(y+dimY/2<tmpDSM.nRows()){
+ //allocate buffer
+ inBuffer.push_back(inBuffer.back());
+ for(int i=0;i<tmpDSM.nCols();++i)
+ inBuffer[inBuffer.size()-1][i]=tmpDSM[y+dimY/2][i];
+ }
+ else{
+ int over=y+dimY/2-tmpDSM.nRows();
+ int index=(inBuffer.size()-1)-over;
+ assert(index>=0);
+ assert(index<inBuffer.size());
+ inBuffer.push_back(inBuffer[index]);
+ }
+ }
+ for(int x=0;x<tmpDSM.nCols();++x){
+ double centerValue=inBuffer[(dimY-1)/2][x];
+ short nmasked=0;
+ std::vector<T> neighbors;
+ for(int j=-(dimY-1)/2;j<=dimY/2;++j){
+ for(int i=-(dimX-1)/2;i<=dimX/2;++i){
+ indexI=x+i;
+ //check if out of bounds
+ if(indexI<0)
+ indexI=-indexI;
+ else if(indexI>=tmpDSM.nCols())
+ indexI=tmpDSM.nCols()-i;
+ if(y+j<0)
+ indexJ=-j;
+ else if(y+j>=tmpDSM.nRows())
+ indexJ=(dimY>2) ? (dimY-1)/2-j : 0;
+ else
+ indexJ=(dimY-1)/2+j;
+ double difference=(centerValue-inBuffer[indexJ][indexI]);
+ if(i||j)//skip centerValue
+ neighbors.push_back(inBuffer[indexJ][indexI]);
+ if(difference>hThreshold)
+ ++nmasked;
+ }
+ }
+ if(nmasked<nlimit){
+ ++nchange;
+ //reset pixel in outputMask
+ outputMask[y][x]=0;
+ }
+ else{
+ //reset pixel height in tmpDSM
+ inBuffer[(dimY-1)/2][x]=stat.mymin(neighbors);
+ }
+ }
+ progress=(1.0+y);
+ progress/=outputMask.nRows();
+ pfnProgress(progress,pszMessage,pProgressArg);
+ }
+ return nchange;
+}
+
template<class T> void Filter2d::shadowDsm(const Vector2d<T>& input, Vector2d<T>& output, double sza, double saa, double pixelSize, short shadowFlag)
{
unsigned int ncols=input.nCols();
diff --git a/src/algorithms/StatFactory.h b/src/algorithms/StatFactory.h
index ed053a7..e6c410b 100644
--- a/src/algorithms/StatFactory.h
+++ b/src/algorithms/StatFactory.h
@@ -754,8 +754,10 @@ template<class T> void StatFactory::percentiles (const std::vector<T>& input, t
inputBin.push_back(*vit);
++vit;
}
- if(inputBin.size())
- output[ibin]=median(inputBin);
+ if(inputBin.size()){
+ /* output[ibin]=median(inputBin); */
+ output[ibin]=mymax(inputBin);
+ }
}
if(!filename.empty()){
std::ofstream outputfile;
diff --git a/src/apps/Makefile.am b/src/apps/Makefile.am
index bdb3053..d8ff05b 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 pkfilterdem pkenhance pkfilterascii pkdsm2shadow pkcomposit pkndvi pkpolygonize pkascii2img pkdiff pksvm pkfssvm pkascii2ogr pkeditogr
+bin_PROGRAMS = pkinfo pkcrop pkreclass pkgetmask pksetmask pkcreatect pkdumpimg pkdumpogr pksieve pkstatascii pkstatogr pkegcs pkextract pkfillnodata pkfilter pkfilterdem pkenhance pkfilterascii pkdsm2shadow pkcomposite pkndvi pkpolygonize pkascii2img pkdiff pksvm pkfssvm pkascii2ogr pkeditogr
# the program to build but not install (the names of the final binaries)
#noinst_PROGRAMS = pkxcorimg pkgeom
@@ -59,7 +59,7 @@ pkenhance_LDADD = $(AM_LDFLAGS) -lgdal
pkfilterascii_SOURCES = pkfilterascii.cc
pkfilterascii_LDADD = $(GSL_LIBS) $(AM_LDFLAGS) -lgsl -lgdal
pkdsm2shadow_SOURCES = pkdsm2shadow.cc
-pkcomposit_SOURCES = pkcomposit.cc
+pkcomposite_SOURCES = pkcomposite.cc
pkndvi_SOURCES = pkndvi.cc
pkpolygonize_SOURCES = pkpolygonize.cc
pkdiff_SOURCES = pkdiff.cc
diff --git a/src/apps/Makefile.in b/src/apps/Makefile.in
index aaead1f..63a5c4c 100644
--- a/src/apps/Makefile.in
+++ b/src/apps/Makefile.in
@@ -57,7 +57,7 @@ bin_PROGRAMS = pkinfo$(EXEEXT) pkcrop$(EXEEXT) pkreclass$(EXEEXT) \
pkstatascii$(EXEEXT) pkstatogr$(EXEEXT) pkegcs$(EXEEXT) \
pkextract$(EXEEXT) pkfillnodata$(EXEEXT) pkfilter$(EXEEXT) \
pkfilterdem$(EXEEXT) pkenhance$(EXEEXT) pkfilterascii$(EXEEXT) \
- pkdsm2shadow$(EXEEXT) pkcomposit$(EXEEXT) pkndvi$(EXEEXT) \
+ pkdsm2shadow$(EXEEXT) pkcomposite$(EXEEXT) pkndvi$(EXEEXT) \
pkpolygonize$(EXEEXT) pkascii2img$(EXEEXT) pkdiff$(EXEEXT) \
pksvm$(EXEEXT) pkfssvm$(EXEEXT) pkascii2ogr$(EXEEXT) \
pkeditogr$(EXEEXT) $(am__EXEEXT_1) $(am__EXEEXT_2) \
@@ -118,10 +118,10 @@ pkascii2ogr_DEPENDENCIES = $(am__DEPENDENCIES_1) $(am__DEPENDENCIES_1) \
$(top_builddir)/src/imageclasses/libimageClasses.la \
$(top_builddir)/src/fileclasses/libfileClasses.la \
$(top_builddir)/src/base/libbase.la
-am_pkcomposit_OBJECTS = pkcomposit.$(OBJEXT)
-pkcomposit_OBJECTS = $(am_pkcomposit_OBJECTS)
-pkcomposit_LDADD = $(LDADD)
-pkcomposit_DEPENDENCIES = $(am__DEPENDENCIES_1) $(am__DEPENDENCIES_1) \
+am_pkcomposite_OBJECTS = pkcomposite.$(OBJEXT)
+pkcomposite_OBJECTS = $(am_pkcomposite_OBJECTS)
+pkcomposite_LDADD = $(LDADD)
+pkcomposite_DEPENDENCIES = $(am__DEPENDENCIES_1) $(am__DEPENDENCIES_1) \
$(top_builddir)/src/algorithms/libalgorithms.la \
$(top_builddir)/src/imageclasses/libimageClasses.la \
$(top_builddir)/src/fileclasses/libfileClasses.la \
@@ -350,7 +350,7 @@ LINK = $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) \
--mode=link $(CCLD) $(AM_CFLAGS) $(CFLAGS) $(AM_LDFLAGS) \
$(LDFLAGS) -o $@
SOURCES = $(pkann_SOURCES) $(pkascii2img_SOURCES) \
- $(pkascii2ogr_SOURCES) $(pkcomposit_SOURCES) \
+ $(pkascii2ogr_SOURCES) $(pkcomposite_SOURCES) \
$(pkcreatect_SOURCES) $(pkcrop_SOURCES) $(pkdiff_SOURCES) \
$(pkdsm2shadow_SOURCES) $(pkdumpimg_SOURCES) \
$(pkdumpogr_SOURCES) $(pkeditogr_SOURCES) $(pkegcs_SOURCES) \
@@ -364,7 +364,7 @@ SOURCES = $(pkann_SOURCES) $(pkascii2img_SOURCES) \
$(pksieve_SOURCES) $(pkstatascii_SOURCES) $(pkstatogr_SOURCES) \
$(pksvm_SOURCES)
DIST_SOURCES = $(am__pkann_SOURCES_DIST) $(pkascii2img_SOURCES) \
- $(pkascii2ogr_SOURCES) $(pkcomposit_SOURCES) \
+ $(pkascii2ogr_SOURCES) $(pkcomposite_SOURCES) \
$(pkcreatect_SOURCES) $(pkcrop_SOURCES) $(pkdiff_SOURCES) \
$(pkdsm2shadow_SOURCES) $(pkdumpimg_SOURCES) \
$(pkdumpogr_SOURCES) $(pkeditogr_SOURCES) $(pkegcs_SOURCES) \
@@ -565,7 +565,7 @@ pkenhance_LDADD = $(AM_LDFLAGS) -lgdal
pkfilterascii_SOURCES = pkfilterascii.cc
pkfilterascii_LDADD = $(GSL_LIBS) $(AM_LDFLAGS) -lgsl -lgdal
pkdsm2shadow_SOURCES = pkdsm2shadow.cc
-pkcomposit_SOURCES = pkcomposit.cc
+pkcomposite_SOURCES = pkcomposite.cc
pkndvi_SOURCES = pkndvi.cc
pkpolygonize_SOURCES = pkpolygonize.cc
pkdiff_SOURCES = pkdiff.cc
@@ -665,9 +665,9 @@ pkascii2img$(EXEEXT): $(pkascii2img_OBJECTS) $(pkascii2img_DEPENDENCIES) $(EXTRA
pkascii2ogr$(EXEEXT): $(pkascii2ogr_OBJECTS) $(pkascii2ogr_DEPENDENCIES) $(EXTRA_pkascii2ogr_DEPENDENCIES)
@rm -f pkascii2ogr$(EXEEXT)
$(CXXLINK) $(pkascii2ogr_OBJECTS) $(pkascii2ogr_LDADD) $(LIBS)
-pkcomposit$(EXEEXT): $(pkcomposit_OBJECTS) $(pkcomposit_DEPENDENCIES) $(EXTRA_pkcomposit_DEPENDENCIES)
- @rm -f pkcomposit$(EXEEXT)
- $(CXXLINK) $(pkcomposit_OBJECTS) $(pkcomposit_LDADD) $(LIBS)
+pkcomposite$(EXEEXT): $(pkcomposite_OBJECTS) $(pkcomposite_DEPENDENCIES) $(EXTRA_pkcomposite_DEPENDENCIES)
+ @rm -f pkcomposite$(EXEEXT)
+ $(CXXLINK) $(pkcomposite_OBJECTS) $(pkcomposite_LDADD) $(LIBS)
pkcreatect$(EXEEXT): $(pkcreatect_OBJECTS) $(pkcreatect_DEPENDENCIES) $(EXTRA_pkcreatect_DEPENDENCIES)
@rm -f pkcreatect$(EXEEXT)
$(CXXLINK) $(pkcreatect_OBJECTS) $(pkcreatect_LDADD) $(LIBS)
@@ -765,7 +765,7 @@ distclean-compile:
@AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/pkann-pkann.Po at am__quote@
@AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/pkascii2img.Po at am__quote@
@AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/pkascii2ogr.Po at am__quote@
- at AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/pkcomposit.Po at am__quote@
+ at AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/pkcomposite.Po at am__quote@
@AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/pkcreatect.Po at am__quote@
@AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/pkcrop.Po at am__quote@
@AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/pkdiff.Po at am__quote@
diff --git a/src/apps/pkann.cc b/src/apps/pkann.cc
index c4eefd4..b199b6d 100644
--- a/src/apps/pkann.cc
+++ b/src/apps/pkann.cc
@@ -1,6 +1,6 @@
/**********************************************************************
-pkclassify_nn.cc: classify raster image using Artificial Neural Network
-Copyright (C) 2008-2012 Pieter Kempeneers
+pkann.cc: classify raster image using Artificial Neural Network
+Copyright (C) 2008-2014 Pieter Kempeneers
This file is part of pktools
@@ -39,13 +39,14 @@ int main(int argc, char *argv[])
//--------------------------- command line options ------------------------------------
Optionpk<string> input_opt("i", "input", "input image");
- Optionpk<string> training_opt("t", "training", "training shape file. A single shape file contains all training features (must be set as: B0, B1, B2,...) for all classes (class numbers identified by label option). Use multiple training files for bootstrap aggregation (alternative to the bag and bsize options, where a random subset is taken from a single training file)");
- Optionpk<string> label_opt("label", "label", "identifier for class label in training shape file.","label");
+ Optionpk<string> training_opt("t", "training", "training vector file. A single vector file contains all training features (must be set as: B0, B1, B2,...) for all classes (class numbers identified by label option). Use multiple training files for bootstrap aggregation (alternative to the bag and bsize options, where a random subset is taken from a single training file)");
+ Optionpk<string> tlayer_opt("tln", "tln", "training layer name(s)");
+ Optionpk<string> label_opt("label", "label", "identifier for class label in training vector file.","label");
Optionpk<unsigned int> balance_opt("bal", "balance", "balance the input data to this number of samples for each class", 0);
Optionpk<bool> random_opt("random", "random", "in case of balance, randomize input data", true);
Optionpk<int> minSize_opt("min", "min", "if number of training pixels is less then min, do not take this class into account (0: consider all classes)", 0);
- Optionpk<double> start_opt("s", "start", "start band sequence number (set to 0)",0);
- Optionpk<double> end_opt("e", "end", "end band sequence number (set to 0 for all bands)", 0);
+ Optionpk<double> start_opt("s", "start", "start band sequence number",0);
+ Optionpk<double> end_opt("e", "end", "end band sequence number (set to 0 to include bands)", 0);
Optionpk<short> band_opt("b", "band", "band index (starting from 0, either use band option or use start to end)");
Optionpk<double> offset_opt("\0", "offset", "offset value for each spectral band input features: refl[band]=(DN[band]-offset[band])/scale[band]", 0.0);
Optionpk<double> scale_opt("\0", "scale", "scale value for each spectral band input features: refl=(DN[band]-offset[band])/scale[band] (use 0 if scale min and max in each band to -1.0 and 1.0)", 0.0);
@@ -83,6 +84,7 @@ int main(int argc, char *argv[])
try{
doProcess=input_opt.retrieveOption(argc,argv);
training_opt.retrieveOption(argc,argv);
+ tlayer_opt.retrieveOption(argc,argv);
label_opt.retrieveOption(argc,argv);
balance_opt.retrieveOption(argc,argv);
random_opt.retrieveOption(argc,argv);
@@ -137,7 +139,7 @@ int main(int argc, char *argv[])
if(mask_opt.size())
cout << "mask filename: " << mask_opt[0] << endl;
if(training_opt.size()){
- cout << "training shape file: " << endl;
+ cout << "training vector file: " << endl;
for(int ifile=0;ifile<training_opt.size();++ifile)
cout << training_opt[ifile] << endl;
}
@@ -207,13 +209,13 @@ int main(int argc, char *argv[])
trainingMap.clear();
trainingPixels.clear();
if(verbose_opt[0]>=1)
- cout << "reading imageShape file " << training_opt[0] << endl;
+ cout << "reading imageVector file " << training_opt[0] << endl;
try{
ImgReaderOgr trainingReaderBag(training_opt[ibag]);
if(band_opt.size())
- totalSamples=trainingReaderBag.readDataImageShape(trainingMap,fields,band_opt,label_opt[0],verbose_opt[0]);
+ totalSamples=trainingReaderBag.readDataImageOgr(trainingMap,fields,band_opt,label_opt[0],tlayer_opt,verbose_opt[0]);
else
- totalSamples=trainingReaderBag.readDataImageShape(trainingMap,fields,start_opt[0],end_opt[0],label_opt[0],verbose_opt[0]);
+ totalSamples=trainingReaderBag.readDataImageOgr(trainingMap,fields,start_opt[0],end_opt[0],label_opt[0],tlayer_opt,verbose_opt[0]);
if(trainingMap.size()<2){
string errorstring="Error: could not read at least two classes from training file, did you provide class labels in training sample (see option label)?";
throw(errorstring);
@@ -576,16 +578,16 @@ int main(int argc, char *argv[])
GDALProgressFunc pfnProgress=GDALTermProgress;
float progress=0;
//-------------------------------- open image file ------------------------------------
- bool refIsRaster=false;
+ bool inputIsRaster=false;
ImgReaderOgr imgReaderOgr;
try{
imgReaderOgr.open(input_opt[0]);
imgReaderOgr.close();
}
catch(string errorString){
- refIsRaster=true;
+ inputIsRaster=true;
}
- if(refIsRaster){
+ if(inputIsRaster){
// if(input_opt[0].find(".shp")==string::npos){
ImgReaderGdal testImage;
try{
@@ -962,9 +964,9 @@ int main(int argc, char *argv[])
classImageBag.close();
classImageOut.close();
}
- else{//classify shape file
+ else{//classify vector file
cm.clearResults();
- //notice that fields have already been set by readDataImageShape (taking into account appropriate bands)
+ //notice that fields have already been set by readDataImageOgr (taking into account appropriate bands)
for(int ivalidation=0;ivalidation<input_opt.size();++ivalidation){
if(output_opt.size())
assert(output_opt.size()==input_opt.size());
@@ -977,18 +979,20 @@ int main(int argc, char *argv[])
if(verbose_opt[0])
std::cout << "opening img writer and copying fields from img reader" << output_opt[ivalidation] << std::endl;
imgWriterOgr.open(output_opt[ivalidation],imgReaderOgr);
- if(verbose_opt[0])
- std::cout << "creating field class" << std::endl;
- if(classValueMap.size())
- imgWriterOgr.createField("class",OFTInteger);
- else
- imgWriterOgr.createField("class",OFTString);
}
if(verbose_opt[0])
cout << "number of layers in input ogr file: " << imgReaderOgr.getLayerCount() << endl;
for(int ilayer=0;ilayer<imgReaderOgr.getLayerCount();++ilayer){
if(verbose_opt[0])
cout << "processing input layer " << ilayer << endl;
+ if(output_opt.size()){
+ if(verbose_opt[0])
+ std::cout << "creating field class" << std::endl;
+ if(classValueMap.size())
+ imgWriterOgr.createField("class",OFTInteger,ilayer);
+ else
+ imgWriterOgr.createField("class",OFTString,ilayer);
+ }
unsigned int nFeatures=imgReaderOgr.getFeatureCount(ilayer);
unsigned int ifeature=0;
progress=0;
@@ -1003,11 +1007,11 @@ int main(int argc, char *argv[])
}
OGRFeature *poDstFeature = NULL;
if(output_opt.size()){
- poDstFeature=imgWriterOgr.createFeature();
+ poDstFeature=imgWriterOgr.createFeature(ilayer);
if( poDstFeature->SetFrom( poFeature, TRUE ) != OGRERR_NONE ){
CPLError( CE_Failure, CPLE_AppDefined,
"Unable to translate feature %d from layer %s.\n",
- poFeature->GetFID(), imgWriterOgr.getLayerName().c_str() );
+ poFeature->GetFID(), imgWriterOgr.getLayerName(ilayer).c_str() );
OGRFeature::DestroyFeature( poFeature );
OGRFeature::DestroyFeature( poDstFeature );
}
@@ -1092,10 +1096,10 @@ int main(int argc, char *argv[])
}
CPLErrorReset();
if(output_opt.size()){
- if(imgWriterOgr.createFeature( poDstFeature ) != OGRERR_NONE){
+ if(imgWriterOgr.createFeature(poDstFeature,ilayer) != OGRERR_NONE){
CPLError( CE_Failure, CPLE_AppDefined,
"Unable to translate feature %d from layer %s.\n",
- poFeature->GetFID(), imgWriterOgr.getLayerName().c_str() );
+ poFeature->GetFID(), imgWriterOgr.getLayerName(ilayer).c_str() );
OGRFeature::DestroyFeature( poDstFeature );
OGRFeature::DestroyFeature( poDstFeature );
}
diff --git a/src/apps/pkascii2img.cc b/src/apps/pkascii2img.cc
index 57bf31c..26a6122 100644
--- a/src/apps/pkascii2img.cc
+++ b/src/apps/pkascii2img.cc
@@ -1,6 +1,6 @@
/**********************************************************************
pkascii2img.cc: program to create raster image based on ascii file
-Copyright (C) 2008-2012 Pieter Kempeneers
+Copyright (C) 2008-2014 Pieter Kempeneers
This file is part of pktools
diff --git a/src/apps/pkascii2ogr.cc b/src/apps/pkascii2ogr.cc
index 2169bca..47989b3 100644
--- a/src/apps/pkascii2ogr.cc
+++ b/src/apps/pkascii2ogr.cc
@@ -1,6 +1,6 @@
/**********************************************************************
pkascii2ogr.cc: program to create vector points or polygons from text file
-Copyright (C) 2008-2012 Pieter Kempeneers
+Copyright (C) 2008-2014 Pieter Kempeneers
This file is part of pktools
diff --git a/src/apps/pkcomposit.cc b/src/apps/pkcomposite.cc
similarity index 97%
rename from src/apps/pkcomposit.cc
rename to src/apps/pkcomposite.cc
index a88a430..7fc6511 100644
--- a/src/apps/pkcomposit.cc
+++ b/src/apps/pkcomposite.cc
@@ -1,6 +1,6 @@
/**********************************************************************
-pkmosaic.cc: program to create mosaic geo-referenced images
-Copyright (C) 2008-2012 Pieter Kempeneers
+pkcomposite.cc: program to mosaic and composite geo-referenced images
+Copyright (C) 2008-2014 Pieter Kempeneers
This file is part of pktools
@@ -51,7 +51,7 @@ int main(int argc, char *argv[])
Optionpk<string> oformat_opt("of", "oformat", "Output image format (see also gdal_translate). Empty string: inherit from input image");
Optionpk<string> colorTable_opt("ct", "ct", "color table (file with 5 columns: id R G B ALFA (0: transparent, 255: solid)");
Optionpk<string> option_opt("co", "co", "Creation option for output file. Multiple options can be specified.");
- Optionpk<unsigned short> resample_opt("r", "resample", "Resampling method (0: nearest neighbour, 1: bi-linear interpolation).", 0);
+ Optionpk<string> resample_opt("r", "resampling-method", "Resampling method (near: nearest neighbour, bilinear: bi-linear interpolation).", "near");
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);
@@ -136,17 +136,19 @@ int main(int argc, char *argv[])
bndnodata_opt.push_back(bndnodata_opt[0]);
}
RESAMPLE theResample;
- switch(resample_opt[0]){
- case(BILINEAR):
- theResample=BILINEAR;
- if(verbose_opt[0])
- cout << "resampling: bilinear interpolation" << endl;
- break;
- default:
+ if(resample_opt[0]=="near"){
theResample=NEAR;
if(verbose_opt[0])
cout << "resampling: nearest neighbour" << endl;
- break;
+ }
+ else if(resample_opt[0]=="bilinear"){
+ theResample=BILINEAR;
+ if(verbose_opt[0])
+ cout << "resampling: bilinear interpolation" << endl;
+ }
+ else{
+ std::cout << "Error: resampling method " << resample_opt[0] << " not supported" << std::endl;
+ exit(1);
}
int nband=0;
@@ -536,7 +538,7 @@ int main(int argc, char *argv[])
double val_current=0;
double val_new=0;
bool readValid=true;
- switch(resample_opt[0]){
+ switch(theResample){
case(BILINEAR):
lowerCol=readCol-0.5;
lowerCol=static_cast<int>(lowerCol);
@@ -606,7 +608,7 @@ int main(int argc, char *argv[])
double ndvi_new=0;
double red_new=0;
double nir_new=0;
- switch(resample_opt[0]){
+ switch(theResample){
case(BILINEAR):
lowerCol=readCol-0.5;
lowerCol=static_cast<int>(lowerCol);
@@ -651,7 +653,7 @@ int main(int argc, char *argv[])
case(crule::minband):
case(crule::validband)://max,min,valid band
val_current=writeBuffer[ruleBand_opt[0]][ib];
- switch(resample_opt[0]){
+ switch(theResample){
case(BILINEAR):
lowerCol=readCol-0.5;
lowerCol=static_cast<int>(lowerCol);
@@ -690,7 +692,7 @@ int main(int argc, char *argv[])
}
break;
case(crule::mode)://max voting (only for Byte images)
- switch(resample_opt[0]){
+ switch(theResample){
case(BILINEAR):
lowerCol=readCol-0.5;
lowerCol=static_cast<int>(lowerCol);
@@ -718,7 +720,7 @@ int main(int argc, char *argv[])
case(crule::mean)://mean value
case(crule::median)://median value
case(crule::sum)://sum value
- switch(resample_opt[0]){
+ switch(theResample){
case(BILINEAR):
lowerCol=readCol-0.5;
lowerCol=static_cast<int>(lowerCol);
@@ -750,7 +752,7 @@ int main(int argc, char *argv[])
break;
case(crule::overwrite):
default:
- switch(resample_opt[0]){
+ switch(theResample){
case(BILINEAR):
lowerCol=readCol-0.5;
lowerCol=static_cast<int>(lowerCol);
@@ -787,7 +789,7 @@ int main(int argc, char *argv[])
case(crule::mean):
case(crule::median):
case(crule::sum):
- switch(resample_opt[0]){
+ switch(theResample){
case(BILINEAR):
lowerCol=readCol-0.5;
lowerCol=static_cast<int>(lowerCol);
@@ -815,7 +817,7 @@ int main(int argc, char *argv[])
++fileBuffer[ib];
break;
case(crule::mode):
- switch(resample_opt[0]){
+ switch(theResample){
case(BILINEAR):
lowerCol=readCol-0.5;
lowerCol=static_cast<int>(lowerCol);
@@ -842,7 +844,7 @@ int main(int argc, char *argv[])
}
break;
default:
- switch(resample_opt[0]){
+ switch(theResample){
case(BILINEAR):
lowerCol=readCol-0.5;
lowerCol=static_cast<int>(lowerCol);
diff --git a/src/apps/pkcreatect.cc b/src/apps/pkcreatect.cc
index 7f14e2b..295cb6f 100644
--- a/src/apps/pkcreatect.cc
+++ b/src/apps/pkcreatect.cc
@@ -1,6 +1,6 @@
/**********************************************************************
pkcreatect.cc: program to create and import colour table to GTiff image
-Copyright (C) 2008-2012 Pieter Kempeneers
+Copyright (C) 2008-2014 Pieter Kempeneers
This file is part of pktools
diff --git a/src/apps/pkcrop.cc b/src/apps/pkcrop.cc
index d27483b..e7466a5 100644
--- a/src/apps/pkcrop.cc
+++ b/src/apps/pkcrop.cc
@@ -1,6 +1,6 @@
/**********************************************************************
pkcrop.cc: perform raster data operations on image such as crop, extract and stack bands
-Copyright (C) 2008-2012 Pieter Kempeneers
+Copyright (C) 2008-2014 Pieter Kempeneers
This file is part of pktools
diff --git a/src/apps/pkdiff.cc b/src/apps/pkdiff.cc
index ae744a2..56415b8 100644
--- a/src/apps/pkdiff.cc
+++ b/src/apps/pkdiff.cc
@@ -1,6 +1,6 @@
/**********************************************************************
pkdiff.cc: program to compare two raster image files
-Copyright (C) 2008-2012 Pieter Kempeneers
+Copyright (C) 2008-2014 Pieter Kempeneers
This file is part of pktools
@@ -262,7 +262,7 @@ int main(int argc, char *argv[])
cout << "created layer" << endl;
cout << "copy fields from " << reference_opt[iref] << endl;
}
- ogrWriter.copyFields(referenceReaderOgr,ilayer);
+ ogrWriter.copyFields(referenceReaderOgr,ilayer,ilayer);
//create extra field for classified label
short theDim=boundary_opt[0];
for(int windowJ=-theDim/2;windowJ<(theDim+1)/2;++windowJ){
diff --git a/src/apps/pkdsm2shadow.cc b/src/apps/pkdsm2shadow.cc
index df62063..ea8d586 100644
--- a/src/apps/pkdsm2shadow.cc
+++ b/src/apps/pkdsm2shadow.cc
@@ -1,6 +1,6 @@
/**********************************************************************
pkdsm2shadow.cc: program to calculate sun shadow based on digital surface model and sun angles
-Copyright (C) 2008-2012 Pieter Kempeneers
+Copyright (C) 2008-2014 Pieter Kempeneers
This file is part of pktools
diff --git a/src/apps/pkdumpimg.cc b/src/apps/pkdumpimg.cc
index 0e40462..3cae28f 100644
--- a/src/apps/pkdumpimg.cc
+++ b/src/apps/pkdumpimg.cc
@@ -1,6 +1,6 @@
/**********************************************************************
pkdumpimg.cc: program to dump image content to ascii or std out
-Copyright (C) 2008-2012 Pieter Kempeneers
+Copyright (C) 2008-2014 Pieter Kempeneers
This file is part of pktools
diff --git a/src/apps/pkdumpogr.cc b/src/apps/pkdumpogr.cc
index 4cf6eff..44b1239 100644
--- a/src/apps/pkdumpogr.cc
+++ b/src/apps/pkdumpogr.cc
@@ -1,6 +1,6 @@
/**********************************************************************
pkdumpogr.cc: dump ogr file to text file or standard output
-Copyright (C) 2008-2012 Pieter Kempeneers
+Copyright (C) 2008-2014 Pieter Kempeneers
This file is part of pktools
diff --git a/src/apps/pkdumpogr.h b/src/apps/pkdumpogr.h
index 4fd4a13..c9a37cd 100644
--- a/src/apps/pkdumpogr.h
+++ b/src/apps/pkdumpogr.h
@@ -1,6 +1,6 @@
/**********************************************************************
pkdumpogr.h: dump ogr file to text file or standard output
-Copyright (C) 2008-2012 Pieter Kempeneers
+Copyright (C) 2008-2014 Pieter Kempeneers
This file is part of pktools
diff --git a/src/apps/pkeditogr.cc b/src/apps/pkeditogr.cc
index b1331be..ba625b8 100644
--- a/src/apps/pkeditogr.cc
+++ b/src/apps/pkeditogr.cc
@@ -1,6 +1,6 @@
/**********************************************************************
pkeditogr.cc: program to edit (rename fields) ogr fil
-Copyright (C) 2008-2013 Pieter Kempeneers
+Copyright (C) 2008-2014 Pieter Kempeneers
This file is part of pktools
diff --git a/src/apps/pkegcs.cc b/src/apps/pkegcs.cc
index 6129a6c..6d772fd 100644
--- a/src/apps/pkegcs.cc
+++ b/src/apps/pkegcs.cc
@@ -1,6 +1,6 @@
/**********************************************************************
pkegcs.cc: Utility for raster files in European Grid Coordinate System
-Copyright (C) 2008-2012 Pieter Kempeneers
+Copyright (C) 2008-2014 Pieter Kempeneers
This file is part of pktools
diff --git a/src/apps/pkenhance.cc b/src/apps/pkenhance.cc
index a518f79..7414156 100644
--- a/src/apps/pkenhance.cc
+++ b/src/apps/pkenhance.cc
@@ -1,6 +1,6 @@
/**********************************************************************
pkenhance.cc: program to enhance raster images: histogram matching
-Copyright (C) 2008-2013 Pieter Kempeneers
+Copyright (C) 2008-2014 Pieter Kempeneers
This file is part of pktools
diff --git a/src/apps/pkextract.cc b/src/apps/pkextract.cc
index 889faf4..2a42071 100644
--- a/src/apps/pkextract.cc
+++ b/src/apps/pkextract.cc
@@ -1,6 +1,6 @@
/**********************************************************************
pkextract.cc: extract pixel values from raster image from a (vector or raster) sample
-Copyright (C) 2008-2012 Pieter Kempeneers
+Copyright (C) 2008-2014 Pieter Kempeneers
This file is part of pktools
@@ -35,18 +35,19 @@ along with pktools. If not, see <http://www.gnu.org/licenses/>.
#endif
namespace rule{
- enum RULE_TYPE {point=0, mean=1, proportion=2, custom=3, minimum=4, maximum=5, maxvote=6, centroid=7, sum=8};
+ enum RULE_TYPE {point=0, mean=1, proportion=2, custom=3, minimum=4, maximum=5, maxvote=6, centroid=7, sum=8, pointOnSurface=9, median=10};
}
using namespace std;
int main(int argc, char *argv[])
{
- Optionpk<string> image_opt("i", "image", "Input image file");
+ Optionpk<string> image_opt("i", "input", "Input image file");
Optionpk<string> sample_opt("s", "sample", "Input sample vector file or class file (e.g. Corine CLC) if class option is set");
+ Optionpk<string> layer_opt("ln", "ln", "layer name(s) in sample (leave empty to select all)");
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. 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. Make sure to set classes if rule is set to maxvote or proportion");
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");
@@ -54,6 +55,7 @@ int main(int argc, char *argv[])
Optionpk<short> geo_opt("g", "geo", "geo coordinates", 1);
Optionpk<short> down_opt("down", "down", "down sampling factor. Can be used to create grid points", 1);
Optionpk<float> threshold_opt("t", "threshold", "threshold for selecting samples (randomly). Provide probability in percentage (>0) or absolute (<0). Use a single threshold for vector sample files. If using raster land cover maps as a sample file, you can provide a threshold value for each class (e.g. -t 80 -t 60). Use value 100 to select all pixels for selected class(es)", 100);
+ Optionpk<float> polythreshold_opt("tp", "thresholdPolygon", "(absolute) threshold for selecting samples in each polygon");
Optionpk<double> min_opt("min", "min", "minimum number of samples to select (0)", 0);
Optionpk<short> boundary_opt("bo", "boundary", "boundary for selecting the sample", 1);
// Optionpk<short> rbox_opt("rb", "rbox", "rectangular boundary box (total width in m) to draw around the selected pixel. Can not combined with class option. Use multiple rbox options for multiple boundary boxes. Use value 0 for no box)", 0);
@@ -65,13 +67,14 @@ int main(int argc, char *argv[])
Optionpk<string> label_opt("cn", "cname", "name of the class label in the output vector file", "label");
Optionpk<bool> polygon_opt("polygon", "polygon", "create OGRPolygon as geometry instead of OGRPoint. Only if sample option is also of polygon type.", false);
Optionpk<int> band_opt("b", "band", "band index to crop. Use -1 to use all bands)", -1);
- Optionpk<string> rule_opt("r", "rule", "rule how to report image information per feature. point (value at each point or at centroid if polygon), centroid, mean (of polygon), proportion, minimum (of polygon), maximum (of polygon), maxvote, sum.", "point");
+ Optionpk<string> rule_opt("r", "rule", "rule how to report image information per feature. point (value at each point or at centroid if polygon), pointOnSurface, centroid, mean (of polygon), median (of polygon), proportion, minimum (of polygon), maximum (of polygon), maxvote, sum.", "point");
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=image_opt.retrieveOption(argc,argv);
sample_opt.retrieveOption(argc,argv);
+ layer_opt.retrieveOption(argc,argv);
mask_opt.retrieveOption(argc,argv);
msknodata_opt.retrieveOption(argc,argv);
class_opt.retrieveOption(argc,argv);
@@ -82,6 +85,7 @@ int main(int argc, char *argv[])
geo_opt.retrieveOption(argc,argv);
down_opt.retrieveOption(argc,argv);
threshold_opt.retrieveOption(argc,argv);
+ polythreshold_opt.retrieveOption(argc,argv);
min_opt.retrieveOption(argc,argv);
boundary_opt.retrieveOption(argc,argv);
// rbox_opt.retrieveOption(argc,argv);
@@ -108,8 +112,10 @@ int main(int argc, char *argv[])
std::map<std::string, rule::RULE_TYPE> ruleMap;
//initialize ruleMap
ruleMap["point"]=rule::point;
+ ruleMap["pointOnSurface"]=rule::pointOnSurface;
ruleMap["centroid"]=rule::centroid;
ruleMap["mean"]=rule::mean;
+ ruleMap["median"]=rule::median;
ruleMap["proportion"]=rule::proportion;
ruleMap["minimum"]=rule::minimum;
ruleMap["maximum"]=rule::maximum;
@@ -248,9 +254,6 @@ int main(int argc, char *argv[])
ImgReaderOgr sampleReaderOgr;
try{
sampleReaderOgr.open(sample_opt[0]);
- //test
- // int nlayer=sampleReaderOgr.getDataSource()->GetLayerCount();
- // cout << "nlayer: " << nlayer << endl;
}
catch(string errorString){
sampleIsRaster=true;
@@ -767,13 +770,19 @@ int main(int argc, char *argv[])
}
//support multiple layers
- int nlayer=sampleReaderOgr.getDataSource()->GetLayerCount();
+ int nlayerRead=sampleReaderOgr.getDataSource()->GetLayerCount();
+ int ilayerWrite=0;
if(verbose_opt[0])
- std::cout << "number of layers: " << nlayer << endl;
+ std::cout << "number of layers: " << nlayerRead << endl;
- for(int ilayer=0;ilayer<nlayer;++ilayer){
+ for(int ilayer=0;ilayer<nlayerRead;++ilayer){
OGRLayer *readLayer=sampleReaderOgr.getLayer(ilayer);
- cout << "processing layer " << readLayer->GetName() << endl;
+ string currentLayername=readLayer->GetName();
+ if(layer_opt.size())
+ if(find(layer_opt.begin(),layer_opt.end(),currentLayername)==layer_opt.end())
+ continue;
+ cout << "processing layer " << currentLayername << endl;
+
readLayer->ResetReading();
OGRLayer *writeLayer;
OGRLayer *writeTestLayer;
@@ -798,46 +807,40 @@ int main(int argc, char *argv[])
}
}
if(verbose_opt[0])
- std::cout << "copy fields" << std::flush << std::endl;
- ogrWriter.copyFields(sampleReaderOgr,ilayer);
+ std::cout << "copy fields from layer " << ilayer << std::flush << std::endl;
+ ogrWriter.copyFields(sampleReaderOgr,ilayer,ilayerWrite);
if(test_opt.size()){
if(verbose_opt[0])
std::cout << "copy fields test writer" << std::flush << std::endl;
- ogrTestWriter.copyFields(sampleReaderOgr,ilayer);
+ ogrTestWriter.copyFields(sampleReaderOgr,ilayer,ilayerWrite);
}
- vector<std::string> fieldnames;
- if(verbose_opt[0])
- std::cout << "get fields" << std::flush << std::endl;
- sampleReaderOgr.getFields(fieldnames);
- assert(fieldnames.size()==ogrWriter.getFieldCount(ilayer));
- map<std::string,double> pointAttributes;
+ // vector<std::string> fieldnames;
+ // if(verbose_opt[0])
+ // std::cout << "get fields" << std::flush << std::endl;
+ // sampleReaderOgr.getFields(fieldnames);
+ // assert(fieldnames.size()==ogrWriter.getFieldCount(ilayerWrite));
+ // map<std::string,double> pointAttributes;
- switch(ruleMap[rule_opt[0]]){
- // switch(rule_opt[0]){
- case(rule::proportion):{//proportion for each class
- assert(class_opt.size());
- for(int iclass=0;iclass<class_opt.size();++iclass){
- ostringstream cs;
- cs << class_opt[iclass];
- ogrWriter.createField(cs.str(),fieldType,ilayer);
+ if(class_opt.size()){
+ switch(ruleMap[rule_opt[0]]){
+ case(rule::proportion):{//proportion for each class
+ for(int iclass=0;iclass<class_opt.size();++iclass){
+ ostringstream cs;
+ cs << class_opt[iclass];
+ ogrWriter.createField(cs.str(),fieldType,ilayerWrite);
+ }
+ break;
}
+ case(rule::custom):
+ case(rule::maxvote):
+ ogrWriter.createField(label_opt[0],fieldType,ilayerWrite);
+ if(test_opt.size())
+ ogrTestWriter.createField(label_opt[0],fieldType,ilayerWrite);
break;
+ }
}
- case(rule::custom):
- case(rule::minimum):
- case(rule::maximum):
- case(rule::maxvote):
- assert(class_opt.size());
- ogrWriter.createField(label_opt[0],fieldType,ilayer);
- if(test_opt.size())
- ogrTestWriter.createField(label_opt[0],fieldType,ilayer);
- break;
- case(rule::point):
- case(rule::mean):
- case(rule::sum):
- case(rule::centroid):
- default:{
+ else{
for(int windowJ=-theDim/2;windowJ<(theDim+1)/2;++windowJ){
for(int windowI=-theDim/2;windowI<(theDim+1)/2;++windowI){
if(disc_opt[0]&&(windowI*windowI+windowJ*windowJ>(theDim/2)*(theDim/2)))
@@ -852,33 +855,16 @@ int main(int argc, char *argv[])
if(verbose_opt[0]>1)
std::cout << "creating field " << fs.str() << std::endl;
- ogrWriter.createField(fs.str(),fieldType,ilayer);
+ ogrWriter.createField(fs.str(),fieldType,ilayerWrite);
if(test_opt.size())
- ogrTestWriter.createField(fs.str(),fieldType,ilayer);
+ ogrTestWriter.createField(fs.str(),fieldType,ilayerWrite);
}
}
}
- break;
- }
}
OGRFeature *readFeature;
unsigned long int ifeature=0;
unsigned long int nfeature=sampleReaderOgr.getFeatureCount();
- // ImgWriterOgr boxWriter;
- // if(rbox_opt[0]>0||cbox_opt[0]>0){
- // assert(bufferOutput_opt.size());
- // assert(test_opt.empty());//not implemented
- // if(verbose_opt[0]>1)
- // std::cout << "opening box writer " << bufferOutput_opt[0] << std::endl;
- // boxWriter.open(bufferOutput_opt[0]);
- // std::string layername="buffer";
- // boxWriter.createLayer(layername, imgReader.getProjection(), wkbPolygon);
- // std::string fieldname="fid";//number of the point
- // if(verbose_opt[0]>1)
- // std::cout << "creating field " << fieldname << std::endl;
- // // ogrWriter.createField(fieldname,OFTInteger,ilayer);
- // boxWriter.createField(fieldname,OFTInteger,ilayer);
- // }
progress=0;
pfnProgress(progress,pszMessage,pProgressArg);
while( (readFeature = readLayer->GetNextFeature()) != NULL ){
@@ -1172,11 +1158,14 @@ int main(int argc, char *argv[])
if(verbose_opt[0]>1)
std::cout << "get centroid point from polygon" << std::endl;
- readPolygon.Centroid(&writeCentroidPoint);
+ if(ruleMap[rule_opt[0]]==rule::pointOnSurface)
+ readPolygon.PointOnSurface(&writeCentroidPoint);
+ else
+ readPolygon.Centroid(&writeCentroidPoint);
double ulx,uly,lrx,lry;
double uli,ulj,lri,lrj;
- if((polygon_opt[0]&&ruleMap[rule_opt[0]]==rule::point)||(ruleMap[rule_opt[0]]==rule::centroid)){
+ if((polygon_opt[0]&&ruleMap[rule_opt[0]]==rule::point)||(ruleMap[rule_opt[0]]==rule::centroid)||(ruleMap[rule_opt[0]]==rule::pointOnSurface)){
ulx=writeCentroidPoint.getX();
uly=writeCentroidPoint.getY();
lrx=ulx;
@@ -1217,263 +1206,223 @@ int main(int argc, char *argv[])
continue;
int nPointPolygon=0;
+
if(polygon_opt[0]){
if(writeTest)
writePolygonFeature = OGRFeature::CreateFeature(writeTestLayer->GetLayerDefn());
else
writePolygonFeature = OGRFeature::CreateFeature(writeLayer->GetLayerDefn());
}
- else if(ruleMap[rule_opt[0]]==rule::mean||ruleMap[rule_opt[0]]==rule::centroid||ruleMap[rule_opt[0]]==rule::sum){
+ else if(ruleMap[rule_opt[0]]!=rule::point){
if(writeTest)
writeCentroidFeature = OGRFeature::CreateFeature(writeTestLayer->GetLayerDefn());
else
writeCentroidFeature = OGRFeature::CreateFeature(writeLayer->GetLayerDefn());
}
- //previously here
- vector<double> polyValues;
- switch(ruleMap[rule_opt[0]]){
- case(rule::point):
- case(rule::mean):
- case(rule::sum):
- case(rule::centroid):
- default:
+ // vector<double> polyValues;
+ Vector2d<double> polyValues;
+ vector<double> polyClassValues;
+
+ if(class_opt.size()){
+ polyClassValues.resize(class_opt.size());
+ //initialize
+ for(int iclass=0;iclass<class_opt.size();++iclass)
+ polyClassValues[iclass]=0;
+ }
+ else
polyValues.resize(nband);
- break;
- case(rule::proportion):
- case(rule::custom):
- case(rule::minimum):
- case(rule::maximum):
- case(rule::maxvote):
- assert(class_opt.size());
- polyValues.resize(class_opt.size());
- break;
+ vector< Vector2d<double> > readValues(nband);
+ for(int iband=0;iband<nband;++iband){
+ int theBand=(band_opt[0]<0)?iband:band_opt[iband];
+ imgReader.readDataBlock(readValues[iband],GDT_Float64,uli,lri,ulj,lrj,theBand);
}
- for(int index=0;index<polyValues.size();++index)
- polyValues[index]=0;
+ //todo: readDataBlock for maskReader...
OGRPoint thePoint;
for(int j=ulj;j<=lrj;++j){
for(int i=uli;i<=lri;++i){
+ //check if within raster image
+ if(i<0||i>=imgReader.nrOfCol())
+ continue;
+ if(j<0||j>=imgReader.nrOfRow())
+ continue;
//check if point is on surface
double x=0;
double y=0;
imgReader.image2geo(i,j,x,y);
thePoint.setX(x);
thePoint.setY(y);
- if(readPolygon.Contains(&thePoint)){
- bool valid=true;
- for(int imask=0;imask<mask_opt.size();++imask){
- double colMask,rowMask;//image coordinates in mask image
- if(mask_opt.size()>1){//multiple masks
- maskReader[imask].geo2image(x,y,colMask,rowMask);
- //nearest neighbour
- rowMask=static_cast<int>(rowMask);
- colMask=static_cast<int>(colMask);
- if(static_cast<int>(colMask)<0||static_cast<int>(colMask)>=maskReader[imask].nrOfCol())
- continue;
- // {
- // cerr << colMask << " out of mask col range!" << std::endl;
- // cerr << x << " " << y << " " << colMask << " " << rowMask << std::endl;
- // assert(static_cast<int>(colMask)>=0&&static_cast<int>(colMask)<maskReader[imask].nrOfCol());
- // }
+
+ if(ruleMap[rule_opt[0]]!=rule::centroid&&!readPolygon.Contains(&thePoint))
+ continue;
+ bool valid=true;
+ for(int imask=0;imask<mask_opt.size();++imask){
+ double colMask,rowMask;//image coordinates in mask image
+ if(mask_opt.size()>1){//multiple masks
+ maskReader[imask].geo2image(x,y,colMask,rowMask);
+ //nearest neighbour
+ rowMask=static_cast<int>(rowMask);
+ colMask=static_cast<int>(colMask);
+ if(static_cast<int>(colMask)<0||static_cast<int>(colMask)>=maskReader[imask].nrOfCol())
+ continue;
- if(static_cast<int>(rowMask)!=static_cast<int>(oldmaskrow[imask])){
- if(static_cast<int>(rowMask)<0||static_cast<int>(rowMask)>=maskReader[imask].nrOfRow())
- continue;
- // {
- // cerr << rowMask << " out of mask row range!" << std::endl;
- // cerr << x << " " << y << " " << colMask << " " << rowMask << std::endl;
- // assert(static_cast<int>(rowMask)>=0&&static_cast<int>(rowMask)<imgReader.nrOfRow());
- // }
- else{
- maskReader[imask].readData(maskBuffer[imask],GDT_Int32,static_cast<int>(rowMask));
- oldmaskrow[imask]=rowMask;
- assert(maskBuffer.size()==maskReader[imask].nrOfBand());
- }
- }
- // char ivalue=0;
- int ivalue=0;
- if(mask_opt.size()==msknodata_opt.size())//one invalid value for each mask
- ivalue=static_cast<int>(msknodata_opt[imask]);
- else//use same invalid value for each mask
- ivalue=static_cast<int>(msknodata_opt[0]);
- if(maskBuffer[imask][colMask]==ivalue){
- valid=false;
- break;
+ if(static_cast<int>(rowMask)!=static_cast<int>(oldmaskrow[imask])){
+ if(static_cast<int>(rowMask)<0||static_cast<int>(rowMask)>=maskReader[imask].nrOfRow())
+ continue;
+ else{
+ maskReader[imask].readData(maskBuffer[imask],GDT_Int32,static_cast<int>(rowMask));
+ oldmaskrow[imask]=rowMask;
+ assert(maskBuffer.size()==maskReader[imask].nrOfBand());
}
}
- else if(maskReader.size()){
- maskReader[0].geo2image(x,y,colMask,rowMask);
- //nearest neighbour
- rowMask=static_cast<int>(rowMask);
- colMask=static_cast<int>(colMask);
- if(static_cast<int>(colMask)<0||static_cast<int>(colMask)>=maskReader[0].nrOfCol())
- continue;
- // {
- // cerr << colMask << " out of mask col range!" << std::endl;
- // cerr << x << " " << y << " " << colMask << " " << rowMask << std::endl;
- // assert(static_cast<int>(colMask)>=0&&static_cast<int>(colMask)<maskReader[0].nrOfCol());
- // }
+ int ivalue=0;
+ if(mask_opt.size()==msknodata_opt.size())//one invalid value for each mask
+ ivalue=static_cast<int>(msknodata_opt[imask]);
+ else//use same invalid value for each mask
+ ivalue=static_cast<int>(msknodata_opt[0]);
+ if(maskBuffer[imask][colMask]==ivalue){
+ valid=false;
+ break;
+ }
+ }
+ else if(maskReader.size()){
+ maskReader[0].geo2image(x,y,colMask,rowMask);
+ //nearest neighbour
+ rowMask=static_cast<int>(rowMask);
+ colMask=static_cast<int>(colMask);
+ if(static_cast<int>(colMask)<0||static_cast<int>(colMask)>=maskReader[0].nrOfCol())
+ continue;
- if(static_cast<int>(rowMask)!=static_cast<int>(oldmaskrow[0])){
- if(static_cast<int>(rowMask)<0||static_cast<int>(rowMask)>=maskReader[0].nrOfRow())
- continue;
- // {
- // cerr << rowMask << " out of mask row range!" << std::endl;
- // cerr << x << " " << y << " " << colMask << " " << rowMask << std::endl;
- // assert(static_cast<int>(rowMask)>=0&&static_cast<int>(rowMask)<imgReader.nrOfRow());
- // }
- else{
- maskReader[0].readData(maskBuffer[0],GDT_Int32,static_cast<int>(rowMask));
- oldmaskrow[0]=rowMask;
- }
+ if(static_cast<int>(rowMask)!=static_cast<int>(oldmaskrow[0])){
+ if(static_cast<int>(rowMask)<0||static_cast<int>(rowMask)>=maskReader[0].nrOfRow())
+ continue;
+ else{
+ maskReader[0].readData(maskBuffer[0],GDT_Int32,static_cast<int>(rowMask));
+ oldmaskrow[0]=rowMask;
}
- for(int ivalue=0;ivalue<msknodata_opt.size();++ivalue){
- if(maskBuffer[0][colMask]==static_cast<int>(msknodata_opt[ivalue])){
- valid=false;
- break;
- }
+ }
+ for(int ivalue=0;ivalue<msknodata_opt.size();++ivalue){
+ if(maskBuffer[0][colMask]==static_cast<int>(msknodata_opt[ivalue])){
+ valid=false;
+ break;
}
}
}
- if(!valid)
- continue;
- else
- validFeature=true;
- //check if within raster image
- if(i<0||i>=imgReader.nrOfCol())
- continue;
- if(j<0||j>=imgReader.nrOfRow())
+ }
+ if(!valid)
+ continue;
+ else
+ validFeature=true;
+ writeRing.addPoint(&thePoint);
+ if(verbose_opt[0]>1)
+ std::cout << "point is on surface:" << thePoint.getX() << "," << thePoint.getY() << std::endl;
+ ++nPointPolygon;
+
+ if(polythreshold_opt.size())
+ if(nPointPolygon>polythreshold_opt[0])
continue;
- writeRing.addPoint(&thePoint);
- if(verbose_opt[0]>1)
- std::cout << "point is on surface:" << thePoint.getX() << "," << thePoint.getY() << std::endl;
- ++nPointPolygon;
- OGRFeature *writePointFeature;
- if(!polygon_opt[0]){
- //create feature
- if(ruleMap[rule_opt[0]]!=rule::mean&&ruleMap[rule_opt[0]]!=rule::centroid&&ruleMap[rule_opt[0]]!=rule::sum){//do not create in case of mean, sum or centroid (only create point at centroid)
- if(writeTest)
- writePointFeature = OGRFeature::CreateFeature(writeTestLayer->GetLayerDefn());
- else
- writePointFeature = OGRFeature::CreateFeature(writeLayer->GetLayerDefn());
- if(verbose_opt[0]>1)
- std::cout << "copying fields from polygons " << sample_opt[0] << std::endl;
- if(writePointFeature->SetFrom(readFeature)!= OGRERR_NONE)
- cerr << "writing feature failed" << std::endl;
- writePointFeature->SetGeometry(&thePoint);
- OGRGeometry *updateGeometry;
- updateGeometry = writePointFeature->GetGeometryRef();
- OGRPoint *poPoint = (OGRPoint *) updateGeometry;
- if(verbose_opt[0]>1)
- std::cout << "write feature has " << writePointFeature->GetFieldCount() << " fields" << std::endl;
- }
+ // throw(nPointPolygon);
+ OGRFeature *writePointFeature;
+ if(!polygon_opt[0]){
+ //create feature
+ if(ruleMap[rule_opt[0]]==rule::point){//do not create in case of mean, median, sum, pointOnSurface or centroid (only create point at centroid)
+ if(writeTest)
+ writePointFeature = OGRFeature::CreateFeature(writeTestLayer->GetLayerDefn());
+ else
+ writePointFeature = OGRFeature::CreateFeature(writeLayer->GetLayerDefn());
+ if(verbose_opt[0]>1)
+ std::cout << "copying fields from polygons " << sample_opt[0] << std::endl;
+ if(writePointFeature->SetFrom(readFeature)!= OGRERR_NONE)
+ cerr << "writing feature failed" << std::endl;
+ writePointFeature->SetGeometry(&thePoint);
+ OGRGeometry *updateGeometry;
+ updateGeometry = writePointFeature->GetGeometryRef();
+ OGRPoint *poPoint = (OGRPoint *) updateGeometry;
+ if(verbose_opt[0]>1)
+ std::cout << "write feature has " << writePointFeature->GetFieldCount() << " fields" << std::endl;
}
- if(verbose_opt[0]>1)
- std::cout << "reading image value within polygon at position " << i << "," << j;
+ }
+ if(class_opt.size()){
+ short value=((readValues[0])[j-ulj])[i-uli];
+ for(int iclass=0;iclass<class_opt.size();++iclass){
+ if(value==class_opt[iclass])
+ polyClassValues[iclass]+=1;
+ }
+ }
+ else{
for(int iband=0;iband<nband;++iband){
int theBand=(band_opt[0]<0)?iband:band_opt[iband];
- double value=0;
- imgReader.readData(value,GDT_Float64,i,j,theBand);
+ assert(j-ulj>=0);
+ assert(j-ulj<readValues[iband].size());
+ assert(i-uli>=0);
+ assert(i-uli<((readValues[iband])[j-ulj]).size());
+ double value=((readValues[iband])[j-ulj])[i-uli];
+ // imgReader.readData(value,GDT_Float64,i,j,theBand);
if(verbose_opt[0]>1)
std::cout << ": " << value << std::endl;
- if(polygon_opt[0]||ruleMap[rule_opt[0]]==rule::mean||ruleMap[rule_opt[0]]==rule::centroid||ruleMap[rule_opt[0]]==rule::sum){
- int iclass=0;
- switch(ruleMap[rule_opt[0]]){
- case(rule::point)://in centroid if polygon_opt==true or all values as points if polygon_opt!=true
- case(rule::centroid):
- default:
- polyValues[iband]=value;
- break;
- case(rule::sum):
- case(rule::mean)://mean as polygon if polygon_opt==true or as point in centroid if polygon_opt!=true
- polyValues[iband]+=value;
- break;
- case(rule::proportion):
- case(rule::custom):
- case(rule::minimum):
- case(rule::maximum):
- case(rule::maxvote):
- for(iclass=0;iclass<class_opt.size();++iclass){
- if(value==class_opt[iclass]){
- assert(polyValues.size()>iclass);
- polyValues[iclass]+=1;//value
- break;
- }
- }
- break;
- }
- }
+ if(polygon_opt[0]||ruleMap[rule_opt[0]]!=rule::point)
+ polyValues[iband].push_back(value);
else{
- // ostringstream fs;
- // if(imgReader.nrOfBand()==1)
- // fs << fieldname_opt[0];
- // else
- // fs << fieldname_opt[0] << theBand;
if(verbose_opt[0]>1)
std::cout << "set field " << fieldname_opt[iband] << " to " << value << std::endl;
switch( fieldType ){
case OFTInteger:
- writePointFeature->SetField(fieldname_opt[iband].c_str(),static_cast<int>(value));
- break;
- case OFTString:
- {
- ostringstream os;
- os << value;
- writePointFeature->SetField(fieldname_opt[iband].c_str(),os.str().c_str());
- break;
- }
case OFTReal:
writePointFeature->SetField(fieldname_opt[iband].c_str(),value);
break;
- case OFTRealList:{
- int fieldIndex=writePointFeature->GetFieldIndex(fieldname_opt[iband].c_str());
- int nCount;
- const double *theList;
- theList=writePointFeature->GetFieldAsDoubleList(fieldIndex,&nCount);
- vector<double> vectorList(nCount+1);
- for(int index=0;index<nCount;++index)
- vectorList[nCount]=value;
- writePointFeature->SetField(fieldIndex,vectorList.size(),&(vectorList[0]));
+ case OFTString:
+ writePointFeature->SetField(fieldname_opt[iband].c_str(),type2string<double>(value).c_str());
break;
- }
+ // case OFTRealList:{
+ // int fieldIndex=writePointFeature->GetFieldIndex(fieldname_opt[iband].c_str());
+ // int nCount;
+ // const double *theList;
+ // theList=writePointFeature->GetFieldAsDoubleList(fieldIndex,&nCount);
+ // vector<double> vectorList(nCount+1);
+ // for(int index=0;index<nCount;++index)
+ // vectorList[nCount]=value;
+ // writePointFeature->SetField(fieldIndex,vectorList.size(),&(vectorList[0]));
+ // break;
+ // }
default://not supported
assert(0);
break;
}
- }
- }
- if(!polygon_opt[0]){
- if(ruleMap[rule_opt[0]]!=rule::mean&&ruleMap[rule_opt[0]]!=rule::centroid&&ruleMap[rule_opt[0]]!=rule::sum){//do not create in case of mean value (only at centroid)
- //write feature
- if(verbose_opt[0]>1)
- std::cout << "creating point feature" << std::endl;
- if(writeTest){
- if(writeTestLayer->CreateFeature( writePointFeature ) != OGRERR_NONE ){
- std::string errorString="Failed to create feature in test shapefile";
- throw(errorString);
- }
+ }//else
+ }//iband
+ }//else (class_opt.size())
+ if(!polygon_opt[0]){
+ if(ruleMap[rule_opt[0]]==rule::point){//do not create in case of mean or median value (only at centroid)
+ //write feature
+ if(verbose_opt[0]>1)
+ std::cout << "creating point feature" << std::endl;
+ if(writeTest){
+ if(writeTestLayer->CreateFeature( writePointFeature ) != OGRERR_NONE ){
+ std::string errorString="Failed to create feature in test shapefile";
+ throw(errorString);
}
- else{
- if(writeLayer->CreateFeature( writePointFeature ) != OGRERR_NONE ){
- std::string errorString="Failed to create feature in shapefile";
- throw(errorString);
- }
+ }
+ else{
+ if(writeLayer->CreateFeature( writePointFeature ) != OGRERR_NONE ){
+ std::string errorString="Failed to create feature in shapefile";
+ throw(errorString);
}
- //destroy feature
- OGRFeature::DestroyFeature( writePointFeature );
- ++ntotalvalid;
- if(verbose_opt[0])
- std::cout << "ntotalvalid(2): " << ntotalvalid << std::endl;
}
+ //destroy feature
+ OGRFeature::DestroyFeature( writePointFeature );
+ ++ntotalvalid;
+ if(verbose_opt[0])
+ std::cout << "ntotalvalid(2): " << ntotalvalid << std::endl;
}
- // ++isample;
}
}
}
- if(polygon_opt[0]||ruleMap[rule_opt[0]]==rule::mean||ruleMap[rule_opt[0]]==rule::centroid||ruleMap[rule_opt[0]]==rule::sum){
+ if(polygon_opt[0]||ruleMap[rule_opt[0]]!=rule::point){
//do not create if no points found within polygon
- if(!nPointPolygon)
+ if(!nPointPolygon){
+ if(verbose_opt[0])
+ cout << "no points found in polygon, continuing" << endl;
continue;
+ }
//add ring to polygon
if(polygon_opt[0]){
writePolygon.addRing(&writeRing);
@@ -1488,7 +1437,7 @@ int main(int argc, char *argv[])
std::cout << "write feature has " << writePolygonFeature->GetFieldCount() << " fields" << std::endl;
//write polygon feature
}
- else{//write mean value of polygon to centroid point (ruleMap[rule_opt[0]]==rule::mean)
+ else{//write value of polygon to centroid point
//create feature
if(verbose_opt[0]>1)
std::cout << "copying fields from polygons " << sample_opt[0] << std::endl;
@@ -1501,258 +1450,202 @@ int main(int argc, char *argv[])
if(verbose_opt[0]>1)
std::cout << "write feature has " << writeCentroidFeature->GetFieldCount() << " fields" << std::endl;
}
- switch(ruleMap[rule_opt[0]]){
- case(rule::point)://value at each point (or at centroid of polygon if line is set)
- default:{
- if(verbose_opt[0])
- std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
- for(int index=0;index<polyValues.size();++index){
- double theValue=polyValues[index];
- // ostringstream fs;
+ if(class_opt.empty()){
+ if(ruleMap[rule_opt[0]]==rule::point){//value at each point (or at centroid of polygon if line is set)
if(verbose_opt[0])
std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
- int theBand=(band_opt[0]<0)?index:band_opt[index];
- // if(nband==1)
- // fs << fieldname_opt[0];
- // else
- // fs << fieldname_opt[0] << theBand;
+ for(int index=0;index<polyValues.size();++index){
+ //test
+ assert(polyValues[index].size()==1);
+ double theValue=polyValues[index].back();
- if(verbose_opt[0]>1)
- std::cout << "set field " << fieldname_opt[index] << " to " << theValue << std::endl;
- switch( fieldType ){
- case OFTInteger:
- if(polygon_opt[0])
- writePolygonFeature->SetField(fieldname_opt[index].c_str(),static_cast<int>(theValue));
- else
- writeCentroidFeature->SetField(fieldname_opt[index].c_str(),static_cast<int>(theValue));
- break;
- case OFTString:
- {
- ostringstream os;
- os << theValue;
+ if(verbose_opt[0])
+ std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
+ int theBand=(band_opt[0]<0)?index:band_opt[index];
+
+ if(verbose_opt[0]>1)
+ std::cout << "set field " << fieldname_opt[index] << " to " << theValue << std::endl;
+ switch( fieldType ){
+ case OFTInteger:
+ case OFTReal:
if(polygon_opt[0])
- writePolygonFeature->SetField(fieldname_opt[index].c_str(),os.str().c_str());
+ writePolygonFeature->SetField(fieldname_opt[index].c_str(),theValue);
else
- writeCentroidFeature->SetField(fieldname_opt[index].c_str(),os.str().c_str());
+ writeCentroidFeature->SetField(fieldname_opt[index].c_str(),theValue);
+ break;
+ case OFTString:
+ if(polygon_opt[0])
+ writePolygonFeature->SetField(fieldname_opt[index].c_str(),type2string<double>(theValue).c_str());
+ else
+ writeCentroidFeature->SetField(fieldname_opt[index].c_str(),type2string<double>(theValue).c_str());
+ break;
+ // case OFTRealList:{
+ // int fieldIndex;
+ // int nCount;
+ // const double *theList;
+ // vector<double> vectorList;
+ // if(polygon_opt[0]){
+ // fieldIndex=writePolygonFeature->GetFieldIndex(fieldname_opt[index].c_str());
+ // theList=writePolygonFeature->GetFieldAsDoubleList(fieldIndex,&nCount);
+ // vectorList.resize(nCount+1);
+ // for(int index=0;index<nCount;++index)
+ // vectorList[index]=theList[index];
+ // vectorList[nCount]=theValue;
+ // writePolygonFeature->SetField(fieldIndex,vectorList.size(),&(vectorList[0]));
+ // }
+ // else{
+ // fieldIndex=writeCentroidFeature->GetFieldIndex(fieldname_opt[index].c_str());
+ // theList=writeCentroidFeature->GetFieldAsDoubleList(fieldIndex,&nCount);
+ // vectorList.resize(nCount+1);
+ // for(int index=0;index<nCount;++index)
+ // vectorList[index]=theList[index];
+ // vectorList[nCount]=theValue;
+ // writeCentroidFeature->SetField(fieldIndex,vectorList.size(),&(vectorList[0]));
+ // }
+ // break;
+ // }
+ default://not supported
+ std::cout << "field type not supported yet..." << std::endl;
break;
}
- case OFTReal:
- if(polygon_opt[0])
- writePolygonFeature->SetField(fieldname_opt[index].c_str(),theValue);
- else
- writeCentroidFeature->SetField(fieldname_opt[index].c_str(),theValue);
- break;
- case OFTRealList:{
- int fieldIndex;
- int nCount;
- const double *theList;
- vector<double> vectorList;
- if(polygon_opt[0]){
- fieldIndex=writePolygonFeature->GetFieldIndex(fieldname_opt[index].c_str());
- theList=writePolygonFeature->GetFieldAsDoubleList(fieldIndex,&nCount);
- vectorList.resize(nCount+1);
- for(int index=0;index<nCount;++index)
- vectorList[index]=theList[index];
- vectorList[nCount]=theValue;
- writePolygonFeature->SetField(fieldIndex,vectorList.size(),&(vectorList[0]));
- }
- else{
- fieldIndex=writeCentroidFeature->GetFieldIndex(fieldname_opt[index].c_str());
- theList=writeCentroidFeature->GetFieldAsDoubleList(fieldIndex,&nCount);
- vectorList.resize(nCount+1);
- for(int index=0;index<nCount;++index)
- vectorList[index]=theList[index];
- vectorList[nCount]=theValue;
- writeCentroidFeature->SetField(fieldIndex,vectorList.size(),&(vectorList[0]));
- }
- break;
- }
- default://not supported
- std::cout << "field type not supported yet..." << std::endl;
- break;
}
}
- break;
- }
- case(rule::mean):
- case(rule::sum):
- case(rule::centroid):{//mean value (written to centroid of polygon if line is not set)
- if(verbose_opt[0])
- std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
- //test
- if(ruleMap[rule_opt[0]]==rule::centroid)
- assert(nPointPolygon<=1);
- for(int index=0;index<polyValues.size();++index){
- double theValue=polyValues[index];
- // ostringstream fs;
- if(ruleMap[rule_opt[0]]==rule::mean)
- theValue/=nPointPolygon;
- int theBand=(band_opt[0]<0)?index:band_opt[index];
- // if(nband==1)
- // fs << fieldname_opt[0];
- // else
- // fs << fieldname_opt[0] << theBand;
- if(verbose_opt[0]>1)
- std::cout << "set field " << fieldname_opt[index] << " to " << theValue << std::endl;
- switch( fieldType ){
- case OFTInteger:
- if(polygon_opt[0])
- writePolygonFeature->SetField(fieldname_opt[index].c_str(),static_cast<int>(theValue));
- else
- writeCentroidFeature->SetField(fieldname_opt[index].c_str(),static_cast<int>(theValue));
- break;
- case OFTString:
- {
- ostringstream os;
- os << theValue;
+ else{//ruleMap[rule_opt[0]] is not rule::point
+ double theValue=0;
+ for(int index=0;index<polyValues.size();++index){
+ if(ruleMap[rule_opt[0]]==rule::mean)
+ theValue=stat.mean(polyValues[index]);
+ else if(ruleMap[rule_opt[0]]==rule::median)
+ theValue=stat.median(polyValues[index]);
+ else if(ruleMap[rule_opt[0]]==rule::sum)
+ theValue=stat.sum(polyValues[index]);
+ else if(ruleMap[rule_opt[0]]==rule::maximum)
+ theValue=stat.mymax(polyValues[index]);
+ else if(ruleMap[rule_opt[0]]==rule::minimum)
+ theValue=stat.mymin(polyValues[index]);
+ else{//rule::pointOnSurface or rule::centroid
+ if(verbose_opt[0])
+ std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
+ assert(nPointPolygon<=1);
+ assert(nPointPolygon==polyValues[index].size());
+ theValue=polyValues[index].back();
+ }
+ if(verbose_opt[0]>1)
+ std::cout << "set field " << fieldname_opt[index] << " to " << theValue << std::endl;
+ switch( fieldType ){
+ case OFTInteger:
+ case OFTReal:
if(polygon_opt[0])
- writePolygonFeature->SetField(fieldname_opt[index].c_str(),os.str().c_str());
+ writePolygonFeature->SetField(fieldname_opt[index].c_str(),theValue);
else
- writeCentroidFeature->SetField(fieldname_opt[index].c_str(),os.str().c_str());
+ writeCentroidFeature->SetField(fieldname_opt[index].c_str(),theValue);
+ break;
+ case OFTString:
+ if(polygon_opt[0])
+ writePolygonFeature->SetField(fieldname_opt[index].c_str(),type2string<double>(theValue).c_str());
+ else
+ writeCentroidFeature->SetField(fieldname_opt[index].c_str(),type2string<double>(theValue).c_str());
+ break;
+ // case OFTRealList:{
+ // int fieldIndex;
+ // int nCount;
+ // const double *theList;
+ // vector<double> vectorList;
+ // if(polygon_opt[0]){
+ // fieldIndex=writePolygonFeature->GetFieldIndex(fieldname_opt[index].c_str());
+ // theList=writePolygonFeature->GetFieldAsDoubleList(fieldIndex,&nCount);
+ // vectorList.resize(nCount+1);
+ // for(int index=0;index<nCount;++index)
+ // vectorList[index]=theList[index];
+ // vectorList[nCount]=theValue;
+ // writePolygonFeature->SetField(fieldIndex,vectorList.size(),&(vectorList[0]));
+ // }
+ // else{
+ // fieldIndex=writeCentroidFeature->GetFieldIndex(fieldname_opt[index].c_str());
+ // theList=writeCentroidFeature->GetFieldAsDoubleList(fieldIndex,&nCount);
+ // vectorList.resize(nCount+1);
+ // for(int index=0;index<nCount;++index)
+ // vectorList[index]=theList[index];
+ // vectorList[nCount]=theValue;
+ // writeCentroidFeature->SetField(fieldIndex,vectorList.size(),&(vectorList[0]));
+ // }
+ // break;
+ //}
+ default://not supported
+ std::cout << "field type not supported yet..." << std::endl;
break;
}
- case OFTReal:
- if(polygon_opt[0])
- writePolygonFeature->SetField(fieldname_opt[index].c_str(),theValue);
- else
- writeCentroidFeature->SetField(fieldname_opt[index].c_str(),theValue);
- break;
- case OFTRealList:{
- int fieldIndex;
- int nCount;
- const double *theList;
- vector<double> vectorList;
- if(polygon_opt[0]){
- fieldIndex=writePolygonFeature->GetFieldIndex(fieldname_opt[index].c_str());
- theList=writePolygonFeature->GetFieldAsDoubleList(fieldIndex,&nCount);
- vectorList.resize(nCount+1);
- for(int index=0;index<nCount;++index)
- vectorList[index]=theList[index];
- vectorList[nCount]=theValue;
- writePolygonFeature->SetField(fieldIndex,vectorList.size(),&(vectorList[0]));
- }
- else{
- fieldIndex=writeCentroidFeature->GetFieldIndex(fieldname_opt[index].c_str());
- theList=writeCentroidFeature->GetFieldAsDoubleList(fieldIndex,&nCount);
- vectorList.resize(nCount+1);
- for(int index=0;index<nCount;++index)
- vectorList[index]=theList[index];
- vectorList[nCount]=theValue;
- writeCentroidFeature->SetField(fieldIndex,vectorList.size(),&(vectorList[0]));
- }
- break;
- }
- default://not supported
- std::cout << "field type not supported yet..." << std::endl;
- break;
}
}
- break;
}
- case(rule::proportion):{//proportion classes
- if(verbose_opt[0])
- std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
- stat.normalize_pct(polyValues);
- // stat.sum(polyValues);
- for(int index=0;index<polyValues.size();++index){
- double theValue=polyValues[index];
- ostringstream fs;
- fs << class_opt[index];
- if(polygon_opt[0])
- writePolygonFeature->SetField(fs.str().c_str(),static_cast<int>(theValue));
- else
- writeCentroidFeature->SetField(fs.str().c_str(),static_cast<int>(theValue));
- }
- break;
- }
- case(rule::custom):{//custom
- assert(polygon_opt[0]);//not implemented for points
- if(verbose_opt[0])
- std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
- stat.normalize_pct(polyValues);
- assert(polyValues.size()==2);//11:broadleaved, 12:coniferous
- if(polyValues[0]>=75)//broadleaved
- writePolygonFeature->SetField(label_opt[0].c_str(),static_cast<int>(11));
- else if(polyValues[1]>=75)//coniferous
- writePolygonFeature->SetField(label_opt[0].c_str(),static_cast<int>(12));
- else if(polyValues[0]>25&&polyValues[1]>25)//mixed
- writePolygonFeature->SetField(label_opt[0].c_str(),static_cast<int>(13));
- else{
- if(verbose_opt[0]){
- std::cout << "No valid value in polyValues..." << std::endl;
- for(int index=0;index<polyValues.size();++index){
- double theValue=polyValues[index];
- std::cout << theValue << " ";
- }
- std::cout << std::endl;
+ else{//class_opt is set
+ if(ruleMap[rule_opt[0]]==rule::proportion){
+ if(verbose_opt[0])
+ std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
+ stat.normalize_pct(polyClassValues);
+ for(int index=0;index<polyValues.size();++index){
+ double theValue=polyClassValues[index];
+ ostringstream fs;
+ fs << class_opt[index];
+ if(polygon_opt[0])
+ writePolygonFeature->SetField(fs.str().c_str(),static_cast<int>(theValue));
+ else
+ writeCentroidFeature->SetField(fs.str().c_str(),static_cast<int>(theValue));
}
- writePolygonFeature->SetField(label_opt[0].c_str(),static_cast<int>(20));
}
- break;
- }
- case(rule::minimum):{
- //minimum of polygon
- assert(polygon_opt[0]);//not implemented for points
- if(verbose_opt[0])
- std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
- //search for min class
- int minClass=stat.mymax(class_opt);
- for(int iclass=0;iclass<class_opt.size();++iclass){
- if(polyValues[iclass]>0){
- if(verbose_opt[0]>1)
- std::cout << class_opt[iclass] << ": " << polyValues[iclass] << std::endl;
- if(class_opt[iclass]<minClass)
- minClass=class_opt[iclass];
+ else if(ruleMap[rule_opt[0]]==rule::custom){
+ assert(polygon_opt[0]);//not implemented for points
+ if(verbose_opt[0])
+ std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
+ stat.normalize_pct(polyClassValues);
+ assert(polyClassValues.size()==2);//11:broadleaved, 12:coniferous
+ if(polyClassValues[0]>=75)//broadleaved
+ writePolygonFeature->SetField(label_opt[0].c_str(),static_cast<int>(11));
+ else if(polyClassValues[1]>=75)//coniferous
+ writePolygonFeature->SetField(label_opt[0].c_str(),static_cast<int>(12));
+ else if(polyClassValues[0]>25&&polyClassValues[1]>25)//mixed
+ writePolygonFeature->SetField(label_opt[0].c_str(),static_cast<int>(13));
+ else{
+ if(verbose_opt[0]){
+ std::cout << "No valid value in polyClassValues..." << std::endl;
+ for(int index=0;index<polyClassValues.size();++index){
+ double theValue=polyClassValues[index];
+ std::cout << theValue << " ";
+ }
+ std::cout << std::endl;
+ }
+ writePolygonFeature->SetField(label_opt[0].c_str(),static_cast<int>(20));
}
}
- if(verbose_opt[0]>0)
- std::cout << "minClass: " << minClass << std::endl;
- writePolygonFeature->SetField(label_opt[0].c_str(),minClass);
- break;
- }
- case(rule::maximum):{
- //maximum of polygon
- assert(polygon_opt[0]);//not implemented for points
- if(verbose_opt[0])
- std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
- //search for max class
- int maxClass=stat.mymin(class_opt);
- for(int iclass=0;iclass<class_opt.size();++iclass){
- if(polyValues[iclass]>0){
- if(verbose_opt[0]>1)
- std::cout << class_opt[iclass] << ": " << polyValues[iclass] << std::endl;
- if(class_opt[iclass]>maxClass)
- maxClass=class_opt[iclass];
- }
+ else if(ruleMap[rule_opt[0]]==rule::maxvote){
+ //maximum votes in polygon
+ if(verbose_opt[0])
+ std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
+ //search for class with maximum votes
+ int maxClass=stat.mymin(class_opt);
+ vector<double>::iterator maxit;
+ maxit=stat.mymax(polyClassValues,polyClassValues.begin(),polyClassValues.end());
+ int maxIndex=distance(polyClassValues.begin(),maxit);
+ maxClass=class_opt[maxIndex];
+ if(verbose_opt[0]>0)
+ std::cout << "maxClass: " << maxClass << std::endl;
+ writePolygonFeature->SetField(label_opt[0].c_str(),maxClass);
}
- if(verbose_opt[0]>0)
- std::cout << "maxClass: " << maxClass << std::endl;
- writePolygonFeature->SetField(label_opt[0].c_str(),maxClass);
- break;
- }
- case(rule::maxvote):{
- //maximum votes in polygon
- assert(polygon_opt[0]);//not implemented for points
- if(verbose_opt[0])
- std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
- //search for class with maximum votes
- int maxClass=stat.mymin(class_opt);
- vector<double>::iterator maxit;
- maxit=stat.mymax(polyValues,polyValues.begin(),polyValues.end());
- int maxIndex=distance(polyValues.begin(),maxit);
- maxClass=class_opt[maxIndex];
- if(verbose_opt[0]>0)
- std::cout << "maxClass: " << maxClass << std::endl;
- writePolygonFeature->SetField(label_opt[0].c_str(),maxClass);
- break;
- }
}
if(polygon_opt[0]){
if(verbose_opt[0]>1)
std::cout << "creating polygon feature" << std::endl;
- if(writeLayer->CreateFeature( writePolygonFeature ) != OGRERR_NONE ){
- std::string errorString="Failed to create polygon feature in shapefile";
- throw(errorString);
+ if(writeTest){
+ if(writeTestLayer->CreateFeature( writePolygonFeature ) != OGRERR_NONE ){
+ std::string errorString="Failed to create polygon feature in shapefile";
+ throw(errorString);
+ }
+ }
+ else{
+ if(writeLayer->CreateFeature( writePolygonFeature ) != OGRERR_NONE ){
+ std::string errorString="Failed to create polygon feature in shapefile";
+ throw(errorString);
+ }
}
OGRFeature::DestroyFeature( writePolygonFeature );
++ntotalvalid;
@@ -1762,9 +1655,19 @@ int main(int argc, char *argv[])
else{
if(verbose_opt[0]>1)
std::cout << "creating point feature in centroid" << std::endl;
- if(writeLayer->CreateFeature( writeCentroidFeature ) != OGRERR_NONE ){
- std::string errorString="Failed to create point feature in shapefile";
- throw(errorString);
+ if(writeTest){
+ if(writeTestLayer->CreateFeature( writeCentroidFeature ) != OGRERR_NONE ){
+ std::string errorString="Failed to create point feature in shapefile";
+ throw(errorString);
+ }
+ }
+ else{
+ //test
+ assert(validFeature);
+ if(writeLayer->CreateFeature( writeCentroidFeature ) != OGRERR_NONE ){
+ std::string errorString="Failed to create point feature in shapefile";
+ throw(errorString);
+ }
}
OGRFeature::DestroyFeature( writeCentroidFeature );
++ntotalvalid;
@@ -1785,6 +1688,7 @@ int main(int argc, char *argv[])
if(verbose_opt[0]>1)
std::cout << "get centroid point from polygon" << std::endl;
+ assert(ruleMap[rule_opt[0]]!=rule::pointOnSurface);//not supported for multipolygons
readPolygon.Centroid(&writeCentroidPoint);
double ulx,uly,lrx,lry;
@@ -1836,45 +1740,49 @@ int main(int argc, char *argv[])
else
writePolygonFeature = OGRFeature::CreateFeature(writeLayer->GetLayerDefn());
}
- else if(ruleMap[rule_opt[0]]==rule::mean||ruleMap[rule_opt[0]]==rule::centroid||ruleMap[rule_opt[0]]==rule::sum){
+ else if(ruleMap[rule_opt[0]]!=rule::point){
if(writeTest)
writeCentroidFeature = OGRFeature::CreateFeature(writeTestLayer->GetLayerDefn());
else
writeCentroidFeature = OGRFeature::CreateFeature(writeLayer->GetLayerDefn());
}
- //previously here
- vector<double> polyValues;
- switch(ruleMap[rule_opt[0]]){
- case(rule::point):
- case(rule::mean):
- case(rule::sum):
- case(rule::centroid):
- default:
+ // vector<double> polyValues;
+ Vector2d<double> polyValues;
+ vector<double> polyClassValues;
+
+ if(class_opt.size()){
+ polyClassValues.resize(class_opt.size());
+ //initialize
+ for(int iclass=0;iclass<class_opt.size();++iclass)
+ polyClassValues[iclass]=0;
+ }
+ else
polyValues.resize(nband);
- break;
- case(rule::proportion):
- case(rule::custom):
- case(rule::minimum):
- case(rule::maximum):
- case(rule::maxvote):
- assert(class_opt.size());
- polyValues.resize(class_opt.size());
- break;
+ vector< Vector2d<double> > readValues(nband);
+ for(int iband=0;iband<nband;++iband){
+ int theBand=(band_opt[0]<0)?iband:band_opt[iband];
+ imgReader.readDataBlock(readValues[iband],GDT_Float64,uli,lri,ulj,lrj,theBand);
}
- for(int index=0;index<polyValues.size();++index)
- polyValues[index]=0;
+ //todo: readDataBlock for maskReader...
OGRPoint thePoint;
for(int j=ulj;j<=lrj;++j){
for(int i=uli;i<=lri;++i){
+ //check if within raster image
+ if(i<0||i>=imgReader.nrOfCol())
+ continue;
+ if(j<0||j>=imgReader.nrOfRow())
+ continue;
//check if point is on surface
double x=0;
double y=0;
imgReader.image2geo(i,j,x,y);
thePoint.setX(x);
thePoint.setY(y);
- if(readPolygon.Contains(&thePoint)){
- bool valid=true;
- for(int imask=0;imask<mask_opt.size();++imask){
+
+ if(ruleMap[rule_opt[0]]!=rule::centroid&&!readPolygon.Contains(&thePoint))
+ continue;
+ bool valid=true;
+ for(int imask=0;imask<mask_opt.size();++imask){
double colMask,rowMask;//image coordinates in mask image
if(mask_opt.size()>1){//multiple masks
maskReader[imask].geo2image(x,y,colMask,rowMask);
@@ -1883,27 +1791,16 @@ int main(int argc, char *argv[])
colMask=static_cast<int>(colMask);
if(static_cast<int>(colMask)<0||static_cast<int>(colMask)>=maskReader[imask].nrOfCol())
continue;
- // {
- // cerr << colMask << " out of mask col range!" << std::endl;
- // cerr << x << " " << y << " " << colMask << " " << rowMask << std::endl;
- // assert(static_cast<int>(colMask)>=0&&static_cast<int>(colMask)<maskReader[imask].nrOfCol());
- // }
if(static_cast<int>(rowMask)!=static_cast<int>(oldmaskrow[imask])){
if(static_cast<int>(rowMask)<0||static_cast<int>(rowMask)>=maskReader[imask].nrOfRow())
continue;
- // {
- // cerr << rowMask << " out of mask row range!" << std::endl;
- // cerr << x << " " << y << " " << colMask << " " << rowMask << std::endl;
- // assert(static_cast<int>(rowMask)>=0&&static_cast<int>(rowMask)<imgReader.nrOfRow());
- // }
else{
maskReader[imask].readData(maskBuffer[imask],GDT_Int32,static_cast<int>(rowMask));
oldmaskrow[imask]=rowMask;
assert(maskBuffer.size()==maskReader[imask].nrOfBand());
}
}
- // char ivalue=0;
int ivalue=0;
if(mask_opt.size()==msknodata_opt.size())//one invalid value for each mask
ivalue=static_cast<int>(msknodata_opt[imask]);
@@ -1921,20 +1818,10 @@ int main(int argc, char *argv[])
colMask=static_cast<int>(colMask);
if(static_cast<int>(colMask)<0||static_cast<int>(colMask)>=maskReader[0].nrOfCol())
continue;
- // {
- // cerr << colMask << " out of mask col range!" << std::endl;
- // cerr << x << " " << y << " " << colMask << " " << rowMask << std::endl;
- // assert(static_cast<int>(colMask)>=0&&static_cast<int>(colMask)<maskReader[0].nrOfCol());
- // }
if(static_cast<int>(rowMask)!=static_cast<int>(oldmaskrow[0])){
if(static_cast<int>(rowMask)<0||static_cast<int>(rowMask)>=maskReader[0].nrOfRow())
continue;
- // {
- // cerr << rowMask << " out of mask row range!" << std::endl;
- // cerr << x << " " << y << " " << colMask << " " << rowMask << std::endl;
- // assert(static_cast<int>(rowMask)>=0&&static_cast<int>(rowMask)<imgReader.nrOfRow());
- // }
else{
maskReader[0].readData(maskBuffer[0],GDT_Int32,static_cast<int>(rowMask));
oldmaskrow[0]=rowMask;
@@ -1952,19 +1839,19 @@ int main(int argc, char *argv[])
continue;
else
validFeature=true;
- //check if within raster image
- if(i<0||i>=imgReader.nrOfCol())
- continue;
- if(j<0||j>=imgReader.nrOfRow())
- continue;
writeRing.addPoint(&thePoint);
if(verbose_opt[0]>1)
std::cout << "point is on surface:" << thePoint.getX() << "," << thePoint.getY() << std::endl;
++nPointPolygon;
+
+ if(polythreshold_opt.size())
+ if(nPointPolygon>polythreshold_opt[0])
+ continue;
+ // throw(nPointPolygon);
OGRFeature *writePointFeature;
if(!polygon_opt[0]){
//create feature
- if(ruleMap[rule_opt[0]]!=rule::mean&&ruleMap[rule_opt[0]]!=rule::centroid&&ruleMap[rule_opt[0]]!=rule::sum){//do not create in case of mean or sum (only create point at centroid)
+ if(ruleMap[rule_opt[0]]==rule::point){//do not create in case of mean, mean or sum (only create point at centroid)
if(writeTest)
writePointFeature = OGRFeature::CreateFeature(writeTestLayer->GetLayerDefn());
else
@@ -1981,82 +1868,57 @@ int main(int argc, char *argv[])
std::cout << "write feature has " << writePointFeature->GetFieldCount() << " fields" << std::endl;
}
}
- if(verbose_opt[0]>1)
- std::cout << "reading image value withinin polygon at position " << i << "," << j;
- for(int iband=0;iband<nband;++iband){
- int theBand=(band_opt[0]<0)?iband:band_opt[iband];
- double value=0;
- imgReader.readData(value,GDT_Float64,i,j,theBand);
- if(verbose_opt[0]>1)
- std::cout << ": " << value << std::endl;
- if(polygon_opt[0]||ruleMap[rule_opt[0]]==rule::mean||ruleMap[rule_opt[0]]==rule::centroid||ruleMap[rule_opt[0]]==rule::sum){
- int iclass=0;
- switch(ruleMap[rule_opt[0]]){
- case(rule::point)://in centroid if polygon_opt==true or all values as points if polygon_opt!=true
- case(rule::centroid):
- default:
- polyValues[iband]=value;
- break;
- case(rule::sum):
- case(rule::mean)://mean or sum polygon if polygon_opt==true or as point in centroid if polygon_opt!=true
- polyValues[iband]+=value;
- break;
- case(rule::proportion):
- case(rule::custom):
- case(rule::minimum):
- case(rule::maximum):
- case(rule::maxvote):
- for(iclass=0;iclass<class_opt.size();++iclass){
- if(value==class_opt[iclass]){
- assert(polyValues.size()>iclass);
- polyValues[iclass]+=1;//value
- break;
- }
- }
- break;
- }
+ if(class_opt.size()){
+ short value=((readValues[0])[j-ulj])[i-uli];
+ for(int iclass=0;iclass<class_opt.size();++iclass){
+ if(value==class_opt[iclass])
+ polyClassValues[iclass]+=1;
}
- else{
- ostringstream fs;
- // if(imgReader.nrOfBand()==1)
- // fs << fieldname_opt[0];
- // else
- // fs << fieldname_opt[0] << theBand;
+ }
+ else{
+ for(int iband=0;iband<nband;++iband){
+ //test
+ assert(j-ulj>=0);
+ assert(j-ulj<readValues[iband].size());
+ assert(i-uli>=0);
+ assert(i-uli<((readValues[iband])[j-ulj]).size());
+ double value=((readValues[iband])[j-ulj])[i-uli];
+ // imgReader.readData(value,GDT_Float64,i,j,theBand);
if(verbose_opt[0]>1)
- std::cout << "set field " << fieldname_opt[iband] << " to " << value << std::endl;
- switch( fieldType ){
- case OFTInteger:
- writePointFeature->SetField(fieldname_opt[iband].c_str(),static_cast<int>(value));
- break;
- case OFTString:
- {
- ostringstream os;
- os << value;
- writePointFeature->SetField(fieldname_opt[iband].c_str(),os.str().c_str());
+ std::cout << ": " << value << std::endl;
+ if(polygon_opt[0]||ruleMap[rule_opt[0]]!=rule::point)
+ polyValues[iband].push_back(value);
+ else{
+ if(verbose_opt[0]>1)
+ std::cout << "set field " << fieldname_opt[iband] << " to " << value << std::endl;
+ switch( fieldType ){
+ case OFTInteger:
+ case OFTReal:
+ writePointFeature->SetField(fieldname_opt[iband].c_str(),value);
+ break;
+ case OFTString:
+ writePointFeature->SetField(fieldname_opt[iband].c_str(),type2string<double>(value).c_str());
+ break;
+ // case OFTRealList:{
+ // int fieldIndex=writePointFeature->GetFieldIndex(fieldname_opt[iband].c_str());
+ // int nCount;
+ // const double *theList;
+ // theList=writePointFeature->GetFieldAsDoubleList(fieldIndex,&nCount);
+ // vector<double> vectorList(nCount+1);
+ // for(int index=0;index<nCount;++index)
+ // vectorList[nCount]=value;
+ // writePointFeature->SetField(fieldIndex,vectorList.size(),&(vectorList[0]));
+ // break;
+ // }
+ default://not supported
+ assert(0);
break;
}
- case OFTReal:
- writePointFeature->SetField(fieldname_opt[iband].c_str(),value);
- break;
- case OFTRealList:{
- int fieldIndex=writePointFeature->GetFieldIndex(fieldname_opt[iband].c_str());
- int nCount;
- const double *theList;
- theList=writePointFeature->GetFieldAsDoubleList(fieldIndex,&nCount);
- vector<double> vectorList(nCount+1);
- for(int index=0;index<nCount;++index)
- vectorList[nCount]=value;
- writePointFeature->SetField(fieldIndex,vectorList.size(),&(vectorList[0]));
- break;
- }
- default://not supported
- assert(0);
- break;
- }
- }
- }
+ }//else
+ }//iband
+ }//else (class_opt.size())
if(!polygon_opt[0]){
- if(ruleMap[rule_opt[0]]!=rule::mean&&ruleMap[rule_opt[0]]!=rule::centroid&&ruleMap[rule_opt[0]]!=rule::sum){//do not create in case of mean value (only at centroid)
+ if(ruleMap[rule_opt[0]]==rule::point){//do not create in case of mean /median value (only at centroid)
//write feature
if(verbose_opt[0]>1)
std::cout << "creating point feature" << std::endl;
@@ -2080,13 +1942,16 @@ int main(int argc, char *argv[])
++ntotalvalid;
if(verbose_opt[0])
std::cout << "ntotalvalid: " << ntotalvalid << std::endl;
- }
}
}
+
//test
if(!validFeature)
continue;
- if(polygon_opt[0]||ruleMap[rule_opt[0]]==rule::mean||ruleMap[rule_opt[0]]==rule::centroid||ruleMap[rule_opt[0]]==rule::sum){
+ if(polygon_opt[0]||ruleMap[rule_opt[0]]!=rule::point){
+ //do not create if no points found within polygon
+ if(!nPointPolygon)
+ continue;
//add ring to polygon
if(polygon_opt[0]){
writePolygon.addRing(&writeRing);
@@ -2101,7 +1966,7 @@ int main(int argc, char *argv[])
std::cout << "write feature has " << writePolygonFeature->GetFieldCount() << " fields" << std::endl;
//write polygon feature
}
- else{//write mean value of polygon to centroid point (ruleMap[rule_opt[0]]==rule::mean)
+ else{//write mean /median value of polygon to centroid point (ruleMap[rule_opt[0]]==rule::mean /median )
//create feature
if(verbose_opt[0]>1)
std::cout << "copying fields from polygons " << sample_opt[0] << std::endl;
@@ -2114,238 +1979,188 @@ int main(int argc, char *argv[])
if(verbose_opt[0]>1)
std::cout << "write feature has " << writeCentroidFeature->GetFieldCount() << " fields" << std::endl;
}
- switch(ruleMap[rule_opt[0]]){
- case(rule::point)://value at each point (or at centroid of polygon if line is set)
- default:{
- if(verbose_opt[0])
- std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
- for(int index=0;index<polyValues.size();++index){
- double theValue=polyValues[index];
- ostringstream fs;
+ if(class_opt.empty()){
+ if(ruleMap[rule_opt[0]]==rule::point){//value at each point (or at centroid of polygon if line is set)
if(verbose_opt[0])
std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
- int theBand=(band_opt[0]<0)?index:band_opt[index];
- // if(nband==1)
- // fs << fieldname_opt[0];
- // else
- // fs << fieldname_opt[0] << theBand;
- if(verbose_opt[0]>1)
- std::cout << "set field " << fieldname_opt[index] << " to " << theValue << std::endl;
- switch( fieldType ){
- case OFTInteger:
- if(polygon_opt[0])
- writePolygonFeature->SetField(fieldname_opt[index].c_str(),static_cast<int>(theValue));
- else
- writeCentroidFeature->SetField(fieldname_opt[index].c_str(),static_cast<int>(theValue));
- break;
- case OFTString:
- {
- ostringstream os;
- os << theValue;
+ for(int index=0;index<polyValues.size();++index){
+ //test
+ assert(polyValues[index].size()==1);
+ double theValue=polyValues[index].back();
+ if(verbose_opt[0])
+ std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
+ int theBand=(band_opt[0]<0)?index:band_opt[index];
+
+ if(verbose_opt[0]>1)
+ std::cout << "set field " << fieldname_opt[index] << " to " << theValue << std::endl;
+ switch( fieldType ){
+ case OFTInteger:
+ case OFTReal:
if(polygon_opt[0])
- writePolygonFeature->SetField(fieldname_opt[index].c_str(),os.str().c_str());
+ writePolygonFeature->SetField(fieldname_opt[index].c_str(),theValue);
else
- writeCentroidFeature->SetField(fieldname_opt[index].c_str(),os.str().c_str());
+ writeCentroidFeature->SetField(fieldname_opt[index].c_str(),theValue);
break;
- }
- case OFTReal:
- if(polygon_opt[0])
- writePolygonFeature->SetField(fieldname_opt[index].c_str(),theValue);
- else
- writeCentroidFeature->SetField(fieldname_opt[index].c_str(),theValue);
- break;
- case OFTRealList:{
- int fieldIndex;
- int nCount;
- const double *theList;
- vector<double> vectorList;
- if(polygon_opt[0]){
- fieldIndex=writePolygonFeature->GetFieldIndex(fieldname_opt[index].c_str());
- theList=writePolygonFeature->GetFieldAsDoubleList(fieldIndex,&nCount);
- vectorList.resize(nCount+1);
- for(int index=0;index<nCount;++index)
- vectorList[index]=theList[index];
- vectorList[nCount]=theValue;
- writePolygonFeature->SetField(fieldIndex,vectorList.size(),&(vectorList[0]));
- }
- else{
- fieldIndex=writeCentroidFeature->GetFieldIndex(fieldname_opt[index].c_str());
- theList=writeCentroidFeature->GetFieldAsDoubleList(fieldIndex,&nCount);
- vectorList.resize(nCount+1);
- for(int index=0;index<nCount;++index)
- vectorList[index]=theList[index];
- vectorList[nCount]=theValue;
- writeCentroidFeature->SetField(fieldIndex,vectorList.size(),&(vectorList[0]));
- }
- break;
- }//case OFTRealList
- }//switch(fieldType)
- }//for index
- break;
- }//case 0 and default
- case(rule::mean):
- case(rule::sum):
- case(rule::centroid):{//mean value (written to centroid of polygon if line is not set)
- if(verbose_opt[0])
- std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
- //test
- if(ruleMap[rule_opt[0]]==rule::centroid)
- assert(nPointPolygon<=1);
- for(int index=0;index<polyValues.size();++index){
- double theValue=polyValues[index];
- // ostringstream fs;
- if(ruleMap[rule_opt[0]]==rule::mean)
- theValue/=nPointPolygon;
- int theBand=(band_opt[0]<0)?index:band_opt[index];
- // if(nband==1)
- // fs << fieldname_opt[0];
- // else
- // fs << fieldname_opt[0] << theBand;
- if(verbose_opt[0]>1)
- std::cout << "set field " << fieldname_opt[index] << " to " << theValue << std::endl;
- switch( fieldType ){
- case OFTInteger:
- if(polygon_opt[0])
- writePolygonFeature->SetField(fieldname_opt[index].c_str(),static_cast<int>(theValue));
- else
- writeCentroidFeature->SetField(fieldname_opt[index].c_str(),static_cast<int>(theValue));
- break;
- case OFTString:
- {
- ostringstream os;
- os << theValue;
+ case OFTString:
if(polygon_opt[0])
- writePolygonFeature->SetField(fieldname_opt[index].c_str(),os.str().c_str());
+ writePolygonFeature->SetField(fieldname_opt[index].c_str(),type2string<double>(theValue).c_str());
else
- writeCentroidFeature->SetField(fieldname_opt[index].c_str(),os.str().c_str());
+ writeCentroidFeature->SetField(fieldname_opt[index].c_str(),type2string<double>(theValue).c_str());
+ break;
+ // case OFTRealList:{
+ // int fieldIndex;
+ // int nCount;
+ // const double *theList;
+ // vector<double> vectorList;
+ // if(polygon_opt[0]){
+ // fieldIndex=writePolygonFeature->GetFieldIndex(fieldname_opt[index].c_str());
+ // theList=writePolygonFeature->GetFieldAsDoubleList(fieldIndex,&nCount);
+ // vectorList.resize(nCount+1);
+ // for(int index=0;index<nCount;++index)
+ // vectorList[index]=theList[index];
+ // vectorList[nCount]=theValue;
+ // writePolygonFeature->SetField(fieldIndex,vectorList.size(),&(vectorList[0]));
+ // }
+ // else{
+ // fieldIndex=writeCentroidFeature->GetFieldIndex(fieldname_opt[index].c_str());
+ // theList=writeCentroidFeature->GetFieldAsDoubleList(fieldIndex,&nCount);
+ // vectorList.resize(nCount+1);
+ // for(int index=0;index<nCount;++index)
+ // vectorList[index]=theList[index];
+ // vectorList[nCount]=theValue;
+ // writeCentroidFeature->SetField(fieldIndex,vectorList.size(),&(vectorList[0]));
+ // }
+ // break;
+ // }
+ default://not supported
+ std::cout << "field type not supported yet..." << std::endl;
break;
}
- case OFTReal:
- if(polygon_opt[0])
- writePolygonFeature->SetField(fieldname_opt[index].c_str(),theValue);
- else
- writeCentroidFeature->SetField(fieldname_opt[index].c_str(),theValue);
- break;
- case OFTRealList:{
- int fieldIndex;
- int nCount;
- const double *theList;
- vector<double> vectorList;
- if(polygon_opt[0]){
- fieldIndex=writePolygonFeature->GetFieldIndex(fieldname_opt[index].c_str());
- theList=writePolygonFeature->GetFieldAsDoubleList(fieldIndex,&nCount);
- vectorList.resize(nCount+1);
- for(int index=0;index<nCount;++index)
- vectorList[index]=theList[index];
- vectorList[nCount]=theValue;
- writePolygonFeature->SetField(fieldIndex,vectorList.size(),&(vectorList[0]));
- }
- else{
- fieldIndex=writeCentroidFeature->GetFieldIndex(fieldname_opt[index].c_str());
- theList=writeCentroidFeature->GetFieldAsDoubleList(fieldIndex,&nCount);
- vectorList.resize(nCount+1);
- for(int index=0;index<nCount;++index)
- vectorList[index]=theList[index];
- vectorList[nCount]=theValue;
- writeCentroidFeature->SetField(fieldIndex,vectorList.size(),&(vectorList[0]));
- }
- break;
}
- }
- }
- break;
- }
- case(rule::proportion):{//proportion classes
- if(verbose_opt[0])
- std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
- stat.normalize_pct(polyValues);
- // stat.sum(polyValues);
- for(int index=0;index<polyValues.size();++index){
- double theValue=polyValues[index];
- ostringstream fs;
- fs << class_opt[index];
- if(polygon_opt[0])
- writePolygonFeature->SetField(fs.str().c_str(),static_cast<int>(theValue));
- else
- writeCentroidFeature->SetField(fs.str().c_str(),static_cast<int>(theValue));
}
- break;
- }
- case(rule::custom):{//custom
- assert(polygon_opt[0]);//not implemented for points
- if(verbose_opt[0])
- std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
- stat.normalize_pct(polyValues);
- assert(polyValues.size()==2);//11:broadleaved, 12:coniferous
- if(polyValues[0]>=75)//broadleaved
- writePolygonFeature->SetField(label_opt[0].c_str(),static_cast<int>(11));
- else if(polyValues[1]>=75)//coniferous
- writePolygonFeature->SetField(label_opt[0].c_str(),static_cast<int>(12));
- else if(polyValues[0]>25&&polyValues[1]>25)//mixed
- writePolygonFeature->SetField(label_opt[0].c_str(),static_cast<int>(13));
- else{
- if(verbose_opt[0]){
- std::cout << "No valid value in polyValues..." << std::endl;
- for(int index=0;index<polyValues.size();++index){
- double theValue=polyValues[index];
- std::cout << theValue << " ";
+ else{//ruleMap[rule_opt[0]] is not rule::point
+ double theValue=0;
+ for(int index=0;index<polyValues.size();++index){
+ if(ruleMap[rule_opt[0]]==rule::mean)
+ theValue=stat.mean(polyValues[index]);
+ else if(ruleMap[rule_opt[0]]==rule::median)
+ theValue=stat.median(polyValues[index]);
+ else if(ruleMap[rule_opt[0]]==rule::sum)
+ theValue=stat.sum(polyValues[index]);
+ else if(ruleMap[rule_opt[0]]==rule::maximum)
+ theValue=stat.mymax(polyValues[index]);
+ else if(ruleMap[rule_opt[0]]==rule::minimum)
+ theValue=stat.mymin(polyValues[index]);
+ else{//rule::pointOnSurface or rule::centroid
+ if(verbose_opt[0])
+ std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
+ assert(nPointPolygon<=1);
+ assert(nPointPolygon==polyValues[index].size());
+ theValue=polyValues[index].back();
+ }
+ if(verbose_opt[0]>1)
+ std::cout << "set field " << fieldname_opt[index] << " to " << theValue << std::endl;
+ switch( fieldType ){
+ case OFTInteger:
+ case OFTReal:
+ if(polygon_opt[0])
+ writePolygonFeature->SetField(fieldname_opt[index].c_str(),theValue);
+ else
+ writeCentroidFeature->SetField(fieldname_opt[index].c_str(),theValue);
+ break;
+ case OFTString:
+ if(polygon_opt[0])
+ writePolygonFeature->SetField(fieldname_opt[index].c_str(),type2string<double>(theValue).c_str());
+ else
+ writeCentroidFeature->SetField(fieldname_opt[index].c_str(),type2string<double>(theValue).c_str());
+ break;
+ // case OFTRealList:{
+ // int fieldIndex;
+ // int nCount;
+ // const double *theList;
+ // vector<double> vectorList;
+ // if(polygon_opt[0]){
+ // fieldIndex=writePolygonFeature->GetFieldIndex(fieldname_opt[index].c_str());
+ // theList=writePolygonFeature->GetFieldAsDoubleList(fieldIndex,&nCount);
+ // vectorList.resize(nCount+1);
+ // for(int index=0;index<nCount;++index)
+ // vectorList[index]=theList[index];
+ // vectorList[nCount]=theValue;
+ // writePolygonFeature->SetField(fieldIndex,vectorList.size(),&(vectorList[0]));
+ // }
+ // else{
+ // fieldIndex=writeCentroidFeature->GetFieldIndex(fieldname_opt[index].c_str());
+ // theList=writeCentroidFeature->GetFieldAsDoubleList(fieldIndex,&nCount);
+ // vectorList.resize(nCount+1);
+ // for(int index=0;index<nCount;++index)
+ // vectorList[index]=theList[index];
+ // vectorList[nCount]=theValue;
+ // writeCentroidFeature->SetField(fieldIndex,vectorList.size(),&(vectorList[0]));
+ // }
+ // break;
+ // }
+ default://not supported
+ std::cout << "field type not supported yet..." << std::endl;
+ break;
}
- std::cout << std::endl;
}
- writePolygonFeature->SetField(label_opt[0].c_str(),static_cast<int>(20));
}
- break;
}
- case(rule::minimum):{//minimum of polygon
- assert(polygon_opt[0]);//not implemented for points
- if(verbose_opt[0])
- std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
- //search for min class
- int minClass=stat.mymax(class_opt);
- for(int iclass=0;iclass<class_opt.size();++iclass){
- if(polyValues[iclass]>0){
- if(verbose_opt[0]>1)
- std::cout << class_opt[iclass] << ": " << polyValues[iclass] << std::endl;
- if(class_opt[iclass]<minClass)
- minClass=class_opt[iclass];
+ else{//class_opt is set
+ if(ruleMap[rule_opt[0]]==rule::proportion){
+ if(verbose_opt[0])
+ std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
+ stat.normalize_pct(polyClassValues);
+ for(int index=0;index<polyValues.size();++index){
+ double theValue=polyClassValues[index];
+ ostringstream fs;
+ fs << class_opt[index];
+ if(polygon_opt[0])
+ writePolygonFeature->SetField(fs.str().c_str(),static_cast<int>(theValue));
+ else
+ writeCentroidFeature->SetField(fs.str().c_str(),static_cast<int>(theValue));
}
}
- if(verbose_opt[0]>0)
- std::cout << "minClass: " << minClass << std::endl;
- writePolygonFeature->SetField(label_opt[0].c_str(),minClass);
- break;
- }
- case(rule::maximum):{//maximum of polygon
- assert(polygon_opt[0]);//not implemented for points
- if(verbose_opt[0])
- std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
- //search for max class
- int maxClass=stat.mymin(class_opt);
- for(int iclass=0;iclass<class_opt.size();++iclass){
- if(polyValues[iclass]>0){
- if(verbose_opt[0]>1)
- std::cout << class_opt[iclass] << ": " << polyValues[iclass] << std::endl;
- if(class_opt[iclass]>maxClass)
- maxClass=class_opt[iclass];
+ else if(ruleMap[rule_opt[0]]==rule::custom){
+ assert(polygon_opt[0]);//not implemented for points
+ if(verbose_opt[0])
+ std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
+ stat.normalize_pct(polyClassValues);
+ assert(polyClassValues.size()==2);//11:broadleaved, 12:coniferous
+ if(polyClassValues[0]>=75)//broadleaved
+ writePolygonFeature->SetField(label_opt[0].c_str(),static_cast<int>(11));
+ else if(polyClassValues[1]>=75)//coniferous
+ writePolygonFeature->SetField(label_opt[0].c_str(),static_cast<int>(12));
+ else if(polyClassValues[0]>25&&polyClassValues[1]>25)//mixed
+ writePolygonFeature->SetField(label_opt[0].c_str(),static_cast<int>(13));
+ else{
+ if(verbose_opt[0]){
+ std::cout << "No valid value in polyClassValues..." << std::endl;
+ for(int index=0;index<polyClassValues.size();++index){
+ double theValue=polyClassValues[index];
+ std::cout << theValue << " ";
+ }
+ std::cout << std::endl;
+ }
+ writePolygonFeature->SetField(label_opt[0].c_str(),static_cast<int>(20));
}
}
- if(verbose_opt[0]>0)
- std::cout << "maxClass: " << maxClass << std::endl;
- writePolygonFeature->SetField(label_opt[0].c_str(),maxClass);
- break;
- }
- case(rule::maxvote):{//maximum votes in polygon
- assert(polygon_opt[0]);//not implemented for points
- if(verbose_opt[0])
- std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
- //search for max votes
- int maxClass=stat.mymin(class_opt);
- vector<double>::iterator maxit;
- maxit=stat.mymax(polyValues,polyValues.begin(),polyValues.end());
- int maxIndex=distance(polyValues.begin(),maxit);
- maxClass=class_opt[maxIndex];
- }
+ else if(ruleMap[rule_opt[0]]==rule::maxvote){
+ //maximum votes in polygon
+ if(verbose_opt[0])
+ std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
+ //search for class with maximum votes
+ int maxClass=stat.mymin(class_opt);
+ vector<double>::iterator maxit;
+ maxit=stat.mymax(polyClassValues,polyClassValues.begin(),polyClassValues.end());
+ int maxIndex=distance(polyClassValues.begin(),maxit);
+ maxClass=class_opt[maxIndex];
+ if(verbose_opt[0]>0)
+ std::cout << "maxClass: " << maxClass << std::endl;
+ writePolygonFeature->SetField(label_opt[0].c_str(),maxClass);
+ }
}
+
if(polygon_opt[0]){
if(verbose_opt[0]>1)
std::cout << "creating polygon feature" << std::endl;
@@ -2405,11 +2220,17 @@ int main(int argc, char *argv[])
std::cout << e << std::endl;
continue;
}
+ catch(int npoint){
+ if(verbose_opt[0])
+ std::cout << "number of points read in polygon: " << npoint << std::endl;
+ continue;
+ }
}//end of getNextFeature
// if(rbox_opt[0]>0||cbox_opt[0]>0)
// boxWriter.close();
progress=1.0;
pfnProgress(progress,pszMessage,pProgressArg);
+ ++ilayerWrite;
}
ogrWriter.close();
if(test_opt.size())
diff --git a/src/apps/pkfillnodata.cc b/src/apps/pkfillnodata.cc
index 6ad1bab..ba6f97b 100644
--- a/src/apps/pkfillnodata.cc
+++ b/src/apps/pkfillnodata.cc
@@ -1,6 +1,6 @@
/**********************************************************************
pkfillnodata.cc: program to fill holes in raster image
-Copyright (C) 2008-2012 Pieter Kempeneers
+Copyright (C) 2008-2014 Pieter Kempeneers
This file is part of pktools
diff --git a/src/apps/pkfilter.cc b/src/apps/pkfilter.cc
index 7d78386..27d38bc 100644
--- a/src/apps/pkfilter.cc
+++ b/src/apps/pkfilter.cc
@@ -1,6 +1,6 @@
/**********************************************************************
pkfilter.cc: program to filter raster images: median, min/max, morphological, filtering
-Copyright (C) 2008-2012 Pieter Kempeneers
+Copyright (C) 2008-2014 Pieter Kempeneers
This file is part of pktools
diff --git a/src/apps/pkfilterascii.cc b/src/apps/pkfilterascii.cc
index d560e4b..f6922a7 100644
--- a/src/apps/pkfilterascii.cc
+++ b/src/apps/pkfilterascii.cc
@@ -1,6 +1,6 @@
/**********************************************************************
pkfilterascii.cc: program to filter data in ASCII file (spectral respons function, dwt)
-Copyright (C) 2008-2013 Pieter Kempeneers
+Copyright (C) 2008-2014 Pieter Kempeneers
This file is part of pktools
diff --git a/src/apps/pkfilterdem.cc b/src/apps/pkfilterdem.cc
index 034d53e..b96e06c 100644
--- a/src/apps/pkfilterdem.cc
+++ b/src/apps/pkfilterdem.cc
@@ -35,8 +35,8 @@ int main(int argc,char **argv) {
Optionpk<std::string> output_opt("o", "output", "Output image file");
Optionpk<std::string> tmpdir_opt("tmp", "tmp", "Temporary directory","/tmp",2);
Optionpk<bool> disc_opt("circ", "circular", "circular disc kernel for dilation and erosion", false);
- Optionpk<string> postFilter_opt("f", "filter", "post processing filter: etew_min, promorph (progressive morphological filter),open,close).");
- Optionpk<double> dim_opt("dim", "dim", "maximum filter kernel size (optionally you can set both initial and maximum filter kernel size", 17);
+ Optionpk<string> postFilter_opt("f", "filter", "post processing filter: vito, etew_min, promorph (progressive morphological filter),open,close).");
+ Optionpk<double> dim_opt("dim", "dim", "maximum filter kernel size (optionally you can set both initial and maximum filter kernel size", 3);
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);
@@ -159,7 +159,131 @@ int main(int argc,char **argv) {
}
unsigned long int nchange=1;
- if(postFilter_opt[0]=="etew_min"){
+ if(postFilter_opt[0]=="vito"){
+ //todo: fill empty pixels
+ // hThreshold_opt.resize(3);
+ // hThreshold_opt[0]=0.7;
+ // hThreshold_opt[1]=0.3;
+ // hThreshold_opt[2]=0.1;
+ vector<int> nlimit(3);
+ nlimit[0]=2;
+ nlimit[1]=3;
+ nlimit[2]=4;
+ nlimit[2]=2;
+ //init finalMask
+ for(int irow=0;irow<tmpData.nRows();++irow)
+ for(int icol=0;icol<tmpData.nCols();++icol)
+ tmpData[irow][icol]=1;
+ for(int iheight=0;iheight<hThreshold_opt.size();++iheight){
+ if(verbose_opt[0])
+ cout << "height: " << hThreshold_opt[iheight] << endl;
+ //todo:replace with binary mask (or short) -> adapt template with T1,T2 in Filter2d
+ Vector2d<double> tmpMask(input.nrOfRow(),input.nrOfCol());
+ for(int irow=0;irow<tmpMask.nRows();++irow)
+ for(int icol=0;icol<tmpMask.nCols();++icol)
+ tmpMask[irow][icol]=1;//1=surface, 0=terrain
+ if(verbose_opt[0])
+ cout << "filtering NWSE" << endl;
+ theFilter.dsm2dtm_nwse(inputData,tmpMask,hThreshold_opt[iheight],nlimit[iheight],dim_opt[0]);
+
+ // //from here
+
+ // Vector2d<double> tmpDSM(inputData);
+ // double noDataValue=0;
+
+ // unsigned long int nchange=0;
+ // int dimX=dim_opt[0];
+ // int dimY=dim_opt[0];
+ // assert(dimX);
+ // assert(dimY);
+ // statfactory::StatFactory stat;
+ // Vector2d<double> inBuffer(dimY,inputData.nCols());
+ // // if(tmpMask.size()!=inputData.nRows())
+ // // tmpMask.resize(inputData.nRows());
+ // int indexI=0;
+ // int indexJ=0;
+ // //initialize last half of inBuffer
+ // for(int j=-(dimY-1)/2;j<=dimY/2;++j){
+ // for(int i=0;i<inputData.nCols();++i)
+ // inBuffer[indexJ][i]=tmpDSM[abs(j)][i];
+ // ++indexJ;
+ // }
+ // for(int y=0;y<tmpDSM.nRows();++y){
+ // if(y){//inBuffer already initialized for y=0
+ // //erase first line from inBuffer
+ // inBuffer.erase(inBuffer.begin());
+ // //read extra line and push back to inBuffer if not out of bounds
+ // if(y+dimY/2<tmpDSM.nRows()){
+ // //allocate buffer
+ // inBuffer.push_back(inBuffer.back());
+ // for(int i=0;i<tmpDSM.nCols();++i)
+ // inBuffer[inBuffer.size()-1][i]=tmpDSM[y+dimY/2][i];
+ // }
+ // else{
+ // int over=y+dimY/2-tmpDSM.nRows();
+ // int index=(inBuffer.size()-1)-over;
+ // assert(index>=0);
+ // assert(index<inBuffer.size());
+ // inBuffer.push_back(inBuffer[index]);
+ // }
+ // }
+ // for(int x=0;x<tmpDSM.nCols();++x){
+ // double centerValue=inBuffer[(dimY-1)/2][x];
+ // short nmasked=0;
+ // std::vector<double> neighbors;
+ // for(int j=-(dimY-1)/2;j<=dimY/2;++j){
+ // for(int i=-(dimX-1)/2;i<=dimX/2;++i){
+ // indexI=x+i;
+ // //check if out of bounds
+ // if(indexI<0)
+ // indexI=-indexI;
+ // else if(indexI>=tmpDSM.nCols())
+ // indexI=tmpDSM.nCols()-i;
+ // if(y+j<0)
+ // indexJ=-j;
+ // else if(y+j>=tmpDSM.nRows())
+ // indexJ=(dimY>2) ? (dimY-1)/2-j : 0;
+ // else
+ // indexJ=(dimY-1)/2+j;
+ // double difference=(centerValue-inBuffer[indexJ][indexI]);
+ // if(i||j)//skip centerValue
+ // neighbors.push_back(inBuffer[indexJ][indexI]);
+ // if(difference>hThreshold_opt[iheight])
+ // ++nmasked;
+ // }
+ // }
+ // if(nmasked<nlimit[iheight]){
+ // ++nchange;
+ // //reset pixel in tmpMask
+ // tmpMask[y][x]=0;
+ // }
+ // else{
+ // //reset pixel height in tmpDSM
+ // inBuffer[(dimY-1)/2][x]=stat.mymin(neighbors);
+ // }
+ // }
+ // }
+ //to here
+
+ tmpData.setMask(tmpMask,0,0);
+ if(verbose_opt[0])
+ cout << "filtering NESW" << endl;
+ theFilter.dsm2dtm_nesw(inputData,tmpMask,hThreshold_opt[iheight],nlimit[iheight],dim_opt[0]);
+ tmpData.setMask(tmpMask,0,0);
+ if(verbose_opt[0])
+ cout << "filtering SENW" << endl;
+ theFilter.dsm2dtm_senw(inputData,tmpMask,hThreshold_opt[iheight],nlimit[iheight],dim_opt[0]);
+ tmpData.setMask(tmpMask,0,0);
+ if(verbose_opt[0])
+ cout << "filtering SWNE" << endl;
+ theFilter.dsm2dtm_swne(inputData,tmpMask,hThreshold_opt[iheight],nlimit[iheight],dim_opt[0]);
+ // set tmpMask to finalMask
+ tmpData.setMask(tmpMask,0,0);
+ }
+ outputData=tmpData;
+ //outputData.setMask(tmpData,1,0);
+ }
+ else 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.
@@ -175,7 +299,7 @@ int main(int argc,char **argv) {
++iteration;
}
}
- else if(postFilter_opt[0]=="promorph"){//todo: instead of number of iterations, define a max dim size
+ else if(postFilter_opt[0]=="promorph"){
//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.
diff --git a/src/apps/pkfsann.cc b/src/apps/pkfsann.cc
index 4a8c939..94fbd3c 100644
--- a/src/apps/pkfsann.cc
+++ b/src/apps/pkfsann.cc
@@ -1,6 +1,6 @@
/**********************************************************************
-pkfs_nn.cc: feature selection for nn classifier
-Copyright (C) 2008-2012 Pieter Kempeneers
+pkfsann.cc: feature selection for nn classifier
+Copyright (C) 2008-2014 Pieter Kempeneers
This file is part of pktools
@@ -197,14 +197,15 @@ int main(int argc, char *argv[])
//--------------------------- command line options ------------------------------------
Optionpk<string> input_opt("i", "input", "input test set (leave empty to perform a cross validation based on training only)");
- Optionpk<string> training_opt("t", "training", "training shape file. A single shape file contains all training features (must be set as: B0, B1, B2,...) for all classes (class numbers identified by label option). Use multiple training files for bootstrap aggregation (alternative to the bag and bsize options, where a random subset is taken from a single training file)");
- Optionpk<string> label_opt("\0", "label", "identifier for class label in training shape file.","label");
+ Optionpk<string> training_opt("t", "training", "training vector file. A single vector file contains all training features (must be set as: B0, B1, B2,...) for all classes (class numbers identified by label option). Use multiple training files for bootstrap aggregation (alternative to the bag and bsize options, where a random subset is taken from a single training file)");
+ Optionpk<string> tlayer_opt("tln", "tln", "training layer name(s)");
+ Optionpk<string> label_opt("\0", "label", "identifier for class label in training vector file.","label");
Optionpk<unsigned short> maxFeatures_opt("n", "nf", "number of features to select (0 to select optimal number, see also ecost option)", 0);
Optionpk<unsigned int> balance_opt("\0", "balance", "balance the input data to this number of samples for each class", 0);
Optionpk<bool> random_opt("random","random", "in case of balance, randomize input data", true);
Optionpk<int> minSize_opt("m", "min", "if number of training pixels is less then min, do not take this class into account", 0);
- Optionpk<double> start_opt("s", "start", "start band sequence number (set to 0)",0);
- Optionpk<double> end_opt("e", "end", "end band sequence number (set to 0 for all bands)", 0);
+ Optionpk<double> start_opt("s", "start", "start band sequence number",0);
+ Optionpk<double> end_opt("e", "end", "end band sequence number (set to 0 to include all bands)", 0);
Optionpk<short> band_opt("b", "band", "band index (starting from 0, either use band option or use start to end)");
Optionpk<double> offset_opt("\0", "offset", "offset value for each spectral band input features: refl[band]=(DN[band]-offset[band])/scale[band]", 0.0);
Optionpk<double> scale_opt("\0", "scale", "scale value for each spectral band input features: refl=(DN[band]-offset[band])/scale[band] (use 0 if scale min and max in each band to -1.0 and 1.0)", 0.0);
@@ -218,6 +219,7 @@ int main(int argc, char *argv[])
doProcess=input_opt.retrieveOption(argc,argv);
training_opt.retrieveOption(argc,argv);
maxFeatures_opt.retrieveOption(argc,argv);
+ tlayer_opt.retrieveOption(argc,argv);
label_opt.retrieveOption(argc,argv);
balance_opt.retrieveOption(argc,argv);
random_opt.retrieveOption(argc,argv);
@@ -261,7 +263,7 @@ int main(int argc, char *argv[])
if(input_opt.size())
cv_opt[0]=0;
if(verbose_opt[0]>=1)
- std::cout << "training shape file: " << training_opt[0] << std::endl;
+ std::cout << "training vector file: " << training_opt[0] << std::endl;
unsigned int totalSamples=0;
unsigned int totalTestSamples=0;
@@ -304,22 +306,22 @@ int main(int argc, char *argv[])
map<string,Vector2d<float> > trainingMap;
map<string,Vector2d<float> > testMap;
if(verbose_opt[0]>=1)
- std::cout << "reading imageShape file " << training_opt[0] << std::endl;
+ std::cout << "reading imageVector file " << training_opt[0] << std::endl;
try{
ImgReaderOgr trainingReader(training_opt[0]);
if(band_opt.size()){
- totalSamples=trainingReader.readDataImageShape(trainingMap,fields,band_opt,label_opt[0],verbose_opt[0]);
+ totalSamples=trainingReader.readDataImageOgr(trainingMap,fields,band_opt,label_opt[0],tlayer_opt,verbose_opt[0]);
if(input_opt.size()){
ImgReaderOgr inputReader(input_opt[0]);
- totalTestSamples=trainingReader.readDataImageShape(testMap,fields,band_opt,label_opt[0],verbose_opt[0]);
+ totalTestSamples=trainingReader.readDataImageOgr(testMap,fields,band_opt,label_opt[0],tlayer_opt,verbose_opt[0]);
inputReader.close();
}
}
else{
- totalSamples=trainingReader.readDataImageShape(trainingMap,fields,start_opt[0],end_opt[0],label_opt[0],verbose_opt[0]);
+ totalSamples=trainingReader.readDataImageOgr(trainingMap,fields,start_opt[0],end_opt[0],label_opt[0],tlayer_opt,verbose_opt[0]);
if(input_opt.size()){
ImgReaderOgr inputReader(input_opt[0]);
- totalTestSamples=trainingReader.readDataImageShape(testMap,fields,start_opt[0],end_opt[0],label_opt[0],verbose_opt[0]);
+ totalTestSamples=trainingReader.readDataImageOgr(testMap,fields,start_opt[0],end_opt[0],label_opt[0],tlayer_opt,verbose_opt[0]);
inputReader.close();
}
}
diff --git a/src/apps/pkfssvm.cc b/src/apps/pkfssvm.cc
index 65c6cab..a314c76 100644
--- a/src/apps/pkfssvm.cc
+++ b/src/apps/pkfssvm.cc
@@ -1,6 +1,6 @@
/**********************************************************************
-pkfs_svm.cc: feature selection for svm classifier
-Copyright (C) 2008-2012 Pieter Kempeneers
+pkfssvm.cc: feature selection for svm classifier
+Copyright (C) 2008-2014 Pieter Kempeneers
This file is part of pktools
@@ -221,14 +221,15 @@ int main(int argc, char *argv[])
//--------------------------- command line options ------------------------------------
Optionpk<string> input_opt("i", "input", "input test set (leave empty to perform a cross validation based on training only)");
- Optionpk<string> training_opt("t", "training", "training shape file. A single shape file contains all training features (must be set as: B0, B1, B2,...) for all classes (class numbers identified by label option).");
- Optionpk<string> label_opt("\0", "label", "identifier for class label in training shape file.","label");
+ Optionpk<string> training_opt("t", "training", "training vector file. A single vector file contains all training features (must be set as: B0, B1, B2,...) for all classes (class numbers identified by label option).");
+ Optionpk<string> tlayer_opt("tln", "tln", "training layer name(s)");
+ Optionpk<string> label_opt("\0", "label", "identifier for class label in training vector file.","label");
Optionpk<unsigned short> maxFeatures_opt("n", "nf", "number of features to select (0 to select optimal number, see also ecost option)", 0);
Optionpk<unsigned int> balance_opt("\0", "balance", "balance the input data to this number of samples for each class", 0);
Optionpk<bool> random_opt("random","random", "in case of balance, randomize input data", true);
Optionpk<int> minSize_opt("m", "min", "if number of training pixels is less then min, do not take this class into account", 0);
- Optionpk<double> start_opt("s", "start", "start band sequence number (set to 0)",0);
- Optionpk<double> end_opt("e", "end", "end band sequence number (set to 0 for all bands)", 0);
+ Optionpk<double> start_opt("s", "start", "start band sequence number",0);
+ Optionpk<double> end_opt("e", "end", "end band sequence number (set to 0 to include all bands)", 0);
Optionpk<short> band_opt("b", "band", "band index (starting from 0, either use band option or use start to end)");
Optionpk<double> offset_opt("\0", "offset", "offset value for each spectral band input features: refl[band]=(DN[band]-offset[band])/scale[band]", 0.0);
Optionpk<double> scale_opt("\0", "scale", "scale value for each spectral band input features: refl=(DN[band]-offset[band])/scale[band] (use 0 if scale min and max in each band to -1.0 and 1.0)", 0.0);
@@ -240,6 +241,7 @@ int main(int argc, char *argv[])
doProcess=input_opt.retrieveOption(argc,argv);
training_opt.retrieveOption(argc,argv);
maxFeatures_opt.retrieveOption(argc,argv);
+ tlayer_opt.retrieveOption(argc,argv);
label_opt.retrieveOption(argc,argv);
balance_opt.retrieveOption(argc,argv);
random_opt.retrieveOption(argc,argv);
@@ -285,7 +287,7 @@ int main(int argc, char *argv[])
if(verbose_opt[0]>=1){
if(input_opt.size())
std::cout << "input filename: " << input_opt[0] << std::endl;
- std::cout << "training shape file: " << std::endl;
+ std::cout << "training vector file: " << std::endl;
for(int ifile=0;ifile<training_opt.size();++ifile)
std::cout << training_opt[ifile] << std::endl;
std::cout << "verbose: " << verbose_opt[0] << std::endl;
@@ -346,18 +348,18 @@ int main(int argc, char *argv[])
try{
ImgReaderOgr trainingReader(training_opt[0]);
if(band_opt.size()){
- totalSamples=trainingReader.readDataImageShape(trainingMap,fields,band_opt,label_opt[0],verbose_opt[0]);
+ totalSamples=trainingReader.readDataImageOgr(trainingMap,fields,band_opt,label_opt[0],tlayer_opt,verbose_opt[0]);
if(input_opt.size()){
ImgReaderOgr inputReader(input_opt[0]);
- totalTestSamples=inputReader.readDataImageShape(testMap,fields,band_opt,label_opt[0],verbose_opt[0]);
+ totalTestSamples=inputReader.readDataImageOgr(testMap,fields,band_opt,label_opt[0],tlayer_opt,verbose_opt[0]);
inputReader.close();
}
}
else{
- totalSamples=trainingReader.readDataImageShape(trainingMap,fields,start_opt[0],end_opt[0],label_opt[0],verbose_opt[0]);
+ totalSamples=trainingReader.readDataImageOgr(trainingMap,fields,start_opt[0],end_opt[0],label_opt[0],tlayer_opt,verbose_opt[0]);
if(input_opt.size()){
ImgReaderOgr inputReader(input_opt[0]);
- totalTestSamples=inputReader.readDataImageShape(testMap,fields,start_opt[0],end_opt[0],label_opt[0],verbose_opt[0]);
+ totalTestSamples=inputReader.readDataImageOgr(testMap,fields,start_opt[0],end_opt[0],label_opt[0],tlayer_opt,verbose_opt[0]);
inputReader.close();
}
}
diff --git a/src/apps/pkgetmask.cc b/src/apps/pkgetmask.cc
index 8c0a861..aa9004e 100644
--- a/src/apps/pkgetmask.cc
+++ b/src/apps/pkgetmask.cc
@@ -1,6 +1,6 @@
/**********************************************************************
pkgetmask.cc: program to create mask image based on values in input raster image
-Copyright (C) 2008-2012 Pieter Kempeneers
+Copyright (C) 2008-2014 Pieter Kempeneers
This file is part of pktools
diff --git a/src/apps/pkinfo.cc b/src/apps/pkinfo.cc
index 6469263..75f1fb1 100644
--- a/src/apps/pkinfo.cc
+++ b/src/apps/pkinfo.cc
@@ -1,6 +1,6 @@
/**********************************************************************
pkinfo.cc: program to retrieve information from raster images
-Copyright (C) 2008-2012 Pieter Kempeneers
+Copyright (C) 2008-2014 Pieter Kempeneers
This file is part of pktools
@@ -31,21 +31,21 @@ int main(int argc, char *argv[])
Optionpk<std::string> input_opt("i","input","Input image file");
Optionpk<bool> bbox_opt("bb", "bbox", "Shows bounding box ", false,0);
Optionpk<bool> bbox_te_opt("te", "te", "Shows bounding box in GDAL format: xmin ymin xmax ymax ", false,0);
- Optionpk<bool> centre_opt("c", "centre", "Image centre in projected X,Y coordinates ", false,0);
- Optionpk<bool> colorTable_opt("ct", "colourtable", "Shows colour table ", false,0);
- Optionpk<bool> samples_opt("ns", "ns", "Number of samples in image ", false,0);
- Optionpk<bool> lines_opt("nl", "nl", "Number of lines in image ", false,0);
+ Optionpk<bool> center_opt("c", "center", "Image center in projected X,Y coordinates ", false,0);
+ Optionpk<bool> colorTable_opt("ct", "colortable", "Shows colour table ", false,0);
+ Optionpk<bool> samples_opt("ns", "nsample", "Number of samples in image ", false,0);
+ Optionpk<bool> lines_opt("nl", "nline", "Number of lines in image ", false,0);
Optionpk<bool> nband_opt("nb", "nband", "Show number of bands in image", false,0);
Optionpk<short> band_opt("b", "band", "Band specific information", 0,0);
Optionpk<bool> dx_opt("dx", "dx", "Gets resolution in x (in m)", false,0);
Optionpk<bool> dy_opt("dy", "dy", "Gets resolution in y (in m)", false,0);
Optionpk<bool> minmax_opt("mm", "minmax", "Shows min and max value of the image ", false,0);
- Optionpk<bool> min_opt("min", "min", "Shows min value of the image ", false,0);
- Optionpk<bool> max_opt("max", "max", "Shows max value of the image ", false,0);
- Optionpk<bool> stat_opt("stats", "stats", "Shows statistics (min,max, mean and stdDev of the image)", false,0);
+ Optionpk<bool> min_opt("min", "minimum", "Shows min value of the image ", false,0);
+ Optionpk<bool> max_opt("max", "maximum", "Shows max value of the image ", false,0);
+ Optionpk<bool> stat_opt("stats", "statistics", "Shows statistics (min,max, mean and stdDev of the image)", false,0);
Optionpk<double> src_min_opt("src_min", "src_min", "Sets minimum for histogram (does not calculate min value: use -mm instead)");
Optionpk<double> src_max_opt("src_max", "src_max", "Sets maximum for histogram (does not calculate min value: use -mm instead)");
- Optionpk<bool> relative_opt("rel", "rel", "Calculates relative histogram in percentage", false,0);
+ Optionpk<bool> relative_opt("rel", "relative", "Calculates relative histogram in percentage", false,0);
Optionpk<bool> projection_opt("a_srs", "a_srs", "Shows projection of the image ", false,0);
Optionpk<bool> geo_opt("geo", "geo", "Gets geotransform ", false,0);
Optionpk<bool> interleave_opt("il", "interleave", "Shows interleave ", false,0);
@@ -54,14 +54,14 @@ int main(int argc, char *argv[])
Optionpk<double> x_opt("x", "xpos", "x pos");
Optionpk<double> y_opt("y", "ypos", "y pos");
Optionpk<bool> read_opt("r", "read", "Reads row y (in projected coordinates if geo option is set, otherwise in image coordinates, 0 based)",false,0);
- Optionpk<bool> refpixel_opt("ref", "ref", "Gets reference pixel (lower left corner of centre of gravity pixel)", false,0);
+ Optionpk<bool> refpixel_opt("ref", "reference", "Gets reference pixel (lower left corner of center of gravity pixel)", false,0);
Optionpk<bool> driver_opt("of", "oformat", "Gets driver description ", false,0);
Optionpk<std::string> extent_opt("e", "extent", "Gets boundary from vector file");
Optionpk<double> ulx_opt("ulx", "ulx", "Upper left x value bounding box");
Optionpk<double> uly_opt("uly", "uly", "Upper left y value bounding box");
Optionpk<double> lrx_opt("lrx", "lrx", "Lower right x value bounding box");
Optionpk<double> lry_opt("lry", "lry", "Lower right y value bounding box");
- Optionpk<bool> hist_opt("hist", "hist", "Calculates histogram. Use --rel for a relative histogram output. ", false,0);
+ Optionpk<bool> hist_opt("hist", "histogram", "Calculates histogram. Use --rel for a relative histogram output. ", false,0);
Optionpk<unsigned int> nbin_opt("nbin", "nbin", "Number of bins used in histogram. Use 0 for all input values as integers");
Optionpk<bool> type_opt("ot", "otype", "Returns data type", false,0);
Optionpk<bool> description_opt("d", "description", "Returns image description", false,0);
@@ -73,7 +73,7 @@ int main(int argc, char *argv[])
doProcess=input_opt.retrieveOption(argc,argv);
bbox_opt.retrieveOption(argc,argv);
bbox_te_opt.retrieveOption(argc,argv);
- centre_opt.retrieveOption(argc,argv);
+ center_opt.retrieveOption(argc,argv);
colorTable_opt.retrieveOption(argc,argv);
samples_opt.retrieveOption(argc,argv);
lines_opt.retrieveOption(argc,argv);
@@ -161,16 +161,16 @@ int main(int argc, char *argv[])
}
if(filename_opt[0])
std::cout << " --input " << input_opt[ifile] << " ";
- if(centre_opt[0]){
+ if(center_opt[0]){
double theX, theY;
- imgReader.getCentrePos(theX,theY);
+ imgReader.getCenterPos(theX,theY);
std::cout << std::setprecision(12) << " -x " << theX << " -y " << theY << " ";
}
if(refpixel_opt[0]){
assert(band_opt[0]<imgReader.nrOfBand());
Egcs egcs;
double refX,refY;
- //get centre of reference (centre of gravity) pixel in image
+ //get center of reference (center of gravity) pixel in image
imgReader.getRefPix(refX,refY,band_opt[0]);
std::cout << std::setprecision(12) << "-x " << refX << " -y " << refY << std::endl;
egcs.setLevel(egcs.res2level(imgReader.getDeltaX()));
@@ -246,12 +246,12 @@ int main(int argc, char *argv[])
}
}
else
- std::cout << "--ct none ";
+ std::cout << "-ct none ";
}
if(samples_opt[0])
- std::cout << "--ns " << imgReader.nrOfCol() << " ";
+ std::cout << "--nsample " << imgReader.nrOfCol() << " ";
if(lines_opt[0])
- std::cout << "--nl " << imgReader.nrOfRow() << " ";
+ std::cout << "--nline " << imgReader.nrOfRow() << " ";
if(nband_opt[0])
std::cout << "--nband " << imgReader.nrOfBand() << " ";
double minValue=0;
@@ -265,7 +265,7 @@ int main(int argc, char *argv[])
GDALRasterBand* rasterBand;
rasterBand=imgReader.getRasterBand(band_opt[0]);
rasterBand->ComputeStatistics(0,&minValue,&maxValue,&meanValue,&stdDev,pfnProgress,pProgressData);
- std::cout << "--min " << minValue << " --max " << maxValue << " --mean " << meanValue << " --stdDev " << stdDev << " ";
+ std::cout << "-min " << minValue << " -max " << maxValue << " --mean " << meanValue << " --stdDev " << stdDev << " ";
}
if(minmax_opt[0]||min_opt[0]||max_opt[0]){
@@ -279,12 +279,12 @@ int main(int argc, char *argv[])
else
imgReader.getMinMax(minValue,maxValue,band_opt[0],true);
if(minmax_opt[0])
- std::cout << "--min " << minValue << " --max " << maxValue << " ";
+ std::cout << "-min " << minValue << " -max " << maxValue << " ";
else{
if(min_opt[0])
- std::cout << "--min " << minValue << " ";
+ std::cout << "-min " << minValue << " ";
if(max_opt[0])
- std::cout << "--max " << maxValue << " ";
+ std::cout << "-max " << maxValue << " ";
}
}
if(relative_opt[0])
@@ -321,11 +321,11 @@ int main(int argc, char *argv[])
// int minCol,minRow;
// if(src_min_opt.size()){
// assert(band_opt[0]<imgReader.nrOfBand());
- // std::cout << "--min " << imgReader.getMin(minCol, minRow,band_opt[0]);
+ // std::cout << "-min " << imgReader.getMin(minCol, minRow,band_opt[0]);
// }
// if(src_max_opt.size()){
// assert(band_opt[0]<imgReader.nrOfBand());
- // std::cout << "--max " << imgReader.getMax(minCol, minRow,band_opt[0]);
+ // std::cout << "-max " << imgReader.getMax(minCol, minRow,band_opt[0]);
// }
// }
if(projection_opt[0]){
@@ -335,7 +335,7 @@ int main(int argc, char *argv[])
std::cout << " -a_srs none" << " ";
}
if(geo_opt[0]&&!read_opt[0]){
- std::cout << " --geo " << std::setprecision(12) << imgReader.getGeoTransform();
+ std::cout << " -geo " << std::setprecision(12) << imgReader.getGeoTransform();
}
if(interleave_opt[0]){
std::cout << " --interleave " << imgReader.getInterleave() << " ";
diff --git a/src/apps/pklas2img.cc b/src/apps/pklas2img.cc
index 293318a..dad94b0 100644
--- a/src/apps/pklas2img.cc
+++ b/src/apps/pklas2img.cc
@@ -1,6 +1,6 @@
/**********************************************************************
pklas2img.cc: program to create (e.g., DEM) raster image from las files
-Copyright (C) 2008-2012 Pieter Kempeneers
+Copyright (C) 2008-2014 Pieter Kempeneers
This file is part of pktools
@@ -32,15 +32,15 @@ int main(int argc,char **argv) {
Optionpk<string> input_opt("i", "input", "Input las file");
Optionpk<short> nodata_opt("nodata", "nodata", "nodata value to put in image if not valid", 0);
Optionpk<string> attribute_opt("n", "name", "names of the attribute to select: intensity, return, nreturn, z", "z");
- 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("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<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("maxit", "maxit", "Maximum number of iterations in post filter", 5);
+ Optionpk<short> nbin_opt("nbin", "nbin", "Number of percentile bins for calculating percentile height value 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> composite_opt("comp", "comp", "composite for multiple points in cell (min, max, median, mean, sum, first, last, profile (percentile height values), number (point density)). Last: overwrite cells with latest point", "last");
+ Optionpk<string> filter_opt("fir", "filter", "filter las points (first,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("dimx", "dimx", "Dimension X of postFilter", 3);
// Optionpk<short> dimy_opt("dimy", "dimy", "Dimension Y of postFilter", 3);
@@ -65,10 +65,10 @@ int main(int argc,char **argv) {
// invalid_opt.retrieveOption(argc,argv);
nodata_opt.retrieveOption(argc,argv);
attribute_opt.retrieveOption(argc,argv);
- disc_opt.retrieveOption(argc,argv);
- maxSlope_opt.retrieveOption(argc,argv);
- hThreshold_opt.retrieveOption(argc,argv);
- maxIter_opt.retrieveOption(argc,argv);
+ // disc_opt.retrieveOption(argc,argv);
+ // maxSlope_opt.retrieveOption(argc,argv);
+ // hThreshold_opt.retrieveOption(argc,argv);
+ // maxIter_opt.retrieveOption(argc,argv);
nbin_opt.retrieveOption(argc,argv);
returns_opt.retrieveOption(argc,argv);
classes_opt.retrieveOption(argc,argv);
@@ -99,7 +99,7 @@ int main(int argc,char **argv) {
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
}
-
+ //todo: is this needed?
GDALAllRegister();
double dfComplete=0.0;
@@ -140,7 +140,7 @@ int main(int argc,char **argv) {
unsigned long int nPoints=0;
unsigned long int ipoint=0;
for(int iinput=0;iinput<input_opt.size();++iinput){
- assert(input_opt[iinput].find(".las")!=string::npos);
+ // assert(input_opt[iinput].find(".las")!=string::npos);
FileReaderLas lasReader;
try{
lasReader.open(input_opt[iinput]);
@@ -149,6 +149,9 @@ int main(int argc,char **argv) {
cout << errorString << endl;
exit(1);
}
+ catch(...){
+ exit(2);
+ }
nPoints=lasReader.getPointCount();
totalPoints+=nPoints;
@@ -230,8 +233,7 @@ int main(int argc,char **argv) {
for(int icol=0;icol<ncol;++icol)
outputData[irow][icol]=0;
-
- std::cout << "Reading " << input_opt.size() << " las files" << std::endl;
+ cout << "Reading " << input_opt.size() << " point cloud files" << endl;
pfnProgress(progress,pszMessage,pProgressArg);
for(int iinput=0;iinput<input_opt.size();++iinput){
FileReaderLas lasReader;
@@ -242,6 +244,12 @@ int main(int argc,char **argv) {
cout << errorString << endl;
exit(1);
}
+ if(verbose_opt[0]){
+ if(lasReader.isCompressed())
+ cout << "Reading compressed point cloud " << input_opt[iinput]<< endl;
+ else
+ cout << "Reading uncompressed point cloud " << input_opt[iinput] << endl;
+ }
//set bounding filter
// lasReader.addBoundsFilter(minULX,maxULY,maxLRX,minLRY);
//set returns filter
@@ -286,10 +294,10 @@ int main(int argc,char **argv) {
continue;
if((filter_opt[0]=="multiple")&&(thePoint.GetNumberOfReturns()<2))
continue;
- if(filter_opt[0]=="last"){
- if(thePoint.GetReturnNumber()!=thePoint.GetNumberOfReturns())
- continue;
- }
+ if((filter_opt[0]=="last")&&(thePoint.GetReturnNumber()!=thePoint.GetNumberOfReturns()))
+ continue;
+ if((filter_opt[0]=="first")&&(thePoint.GetReturnNumber()!=1))
+ continue;
double dcol,drow;
outputWriter.geo2image(thePoint.GetX(),thePoint.GetY(),dcol,drow);
int icol=static_cast<int>(dcol);
diff --git a/src/apps/pkndvi.cc b/src/apps/pkndvi.cc
index 4ac2b1e..210afd5 100644
--- a/src/apps/pkndvi.cc
+++ b/src/apps/pkndvi.cc
@@ -1,6 +1,6 @@
/**********************************************************************
pkndvi.cc: program to calculate vegetation index image
-Copyright (C) 2008-2012 Pieter Kempeneers
+Copyright (C) 2008-2014 Pieter Kempeneers
This file is part of pktools
diff --git a/src/apps/pkoptsvm.cc b/src/apps/pkoptsvm.cc
index b1422bc..8980d0a 100644
--- a/src/apps/pkoptsvm.cc
+++ b/src/apps/pkoptsvm.cc
@@ -1,6 +1,6 @@
/**********************************************************************
-pkopt_svm.cc: program to optimize parameters for SVM classification
-Copyright (C) 2008-2013 Pieter Kempeneers
+pkoptsvm.cc: program to optimize parameters for SVM classification
+Copyright (C) 2008-2014 Pieter Kempeneers
This file is part of pktools
@@ -231,14 +231,15 @@ int main(int argc, char *argv[])
map<short,int> reclassMap;
vector<int> vreclass;
Optionpk<string> input_opt("i", "input", "input image");
- Optionpk<string> training_opt("t", "training", "training shape file. A single shape file contains all training features (must be set as: B0, B1, B2,...) for all classes (class numbers identified by label option).");
- Optionpk<string> label_opt("\0", "label", "identifier for class label in training shape file.","label");
+ Optionpk<string> training_opt("t", "training", "training vector file. A single vector file contains all training features (must be set as: B0, B1, B2,...) for all classes (class numbers identified by label option).");
+ Optionpk<string> tlayer_opt("tln", "tln", "training layer name(s)");
+ Optionpk<string> label_opt("\0", "label", "identifier for class label in training vector file.","label");
// Optionpk<unsigned short> reclass_opt("\0", "rc", "reclass code (e.g. --rc=12 --rc=23 to reclass first two classes to 12 and 23 resp.).", 0);
Optionpk<unsigned int> balance_opt("\0", "balance", "balance the input data to this number of samples for each class", 0);
Optionpk<bool> random_opt("random","random", "in case of balance, randomize input data", true);
Optionpk<int> minSize_opt("m", "min", "if number of training pixels is less then min, do not take this class into account", 0);
- Optionpk<double> start_opt("s", "start", "start band sequence number (set to 0)",0);
- Optionpk<double> end_opt("e", "end", "end band sequence number (set to 0 for all bands)", 0);
+ Optionpk<double> start_opt("s", "start", "start band sequence number",0);
+ Optionpk<double> end_opt("e", "end", "end band sequence number (set to 0 to include all bands)", 0);
Optionpk<short> band_opt("b", "band", "band index (starting from 0, either use band option or use start to end)");
Optionpk<double> offset_opt("\0", "offset", "offset value for each spectral band input features: refl[band]=(DN[band]-offset[band])/scale[band]", 0.0);
Optionpk<double> scale_opt("\0", "scale", "scale value for each spectral band input features: refl=(DN[band]-offset[band])/scale[band] (use 0 if scale min and max in each band to -1.0 and 1.0)", 0.0);
@@ -253,6 +254,7 @@ int main(int argc, char *argv[])
try{
doProcess=input_opt.retrieveOption(argc,argv);
training_opt.retrieveOption(argc,argv);
+ tlayer_opt.retrieveOption(argc,argv);
label_opt.retrieveOption(argc,argv);
// reclass_opt.retrieveOption(argc,argv);
balance_opt.retrieveOption(argc,argv);
@@ -301,7 +303,7 @@ int main(int argc, char *argv[])
if(verbose_opt[0]>=1){
if(input_opt.size())
std::cout << "input filename: " << input_opt[0] << std::endl;
- std::cout << "training shape file: " << std::endl;
+ std::cout << "training vector file: " << std::endl;
for(int ifile=0;ifile<training_opt.size();++ifile)
std::cout << training_opt[ifile] << std::endl;
std::cout << "verbose: " << verbose_opt[0] << std::endl;
@@ -355,18 +357,18 @@ int main(int argc, char *argv[])
try{
ImgReaderOgr trainingReader(training_opt[0]);
if(band_opt.size()){
- totalSamples=trainingReader.readDataImageShape(trainingMap,fields,band_opt,label_opt[0],verbose_opt[0]);
+ totalSamples=trainingReader.readDataImageOgr(trainingMap,fields,band_opt,label_opt[0],tlayer_opt,verbose_opt[0]);
if(input_opt.size()){
ImgReaderOgr inputReader(input_opt[0]);
- totalTestSamples=inputReader.readDataImageShape(testMap,fields,band_opt,label_opt[0],verbose_opt[0]);
+ totalTestSamples=inputReader.readDataImageOgr(testMap,fields,band_opt,label_opt[0],tlayer_opt,verbose_opt[0]);
inputReader.close();
}
}
else{
- totalSamples=trainingReader.readDataImageShape(trainingMap,fields,start_opt[0],end_opt[0],label_opt[0],verbose_opt[0]);
+ totalSamples=trainingReader.readDataImageOgr(trainingMap,fields,start_opt[0],end_opt[0],label_opt[0],tlayer_opt,verbose_opt[0]);
if(input_opt.size()){
ImgReaderOgr inputReader(input_opt[0]);
- totalTestSamples=inputReader.readDataImageShape(testMap,fields,start_opt[0],end_opt[0],label_opt[0],verbose_opt[0]);
+ totalTestSamples=inputReader.readDataImageOgr(testMap,fields,start_opt[0],end_opt[0],label_opt[0],tlayer_opt,verbose_opt[0]);
inputReader.close();
}
trainingReader.close();
diff --git a/src/apps/pkpolygonize.cc b/src/apps/pkpolygonize.cc
index 55b0b5c..ef2e85b 100644
--- a/src/apps/pkpolygonize.cc
+++ b/src/apps/pkpolygonize.cc
@@ -1,6 +1,6 @@
/**********************************************************************
pkpolygonize.cc: program to make vector file from raster image
-Copyright (C) 2008-2012 Pieter Kempeneers
+Copyright (C) 2008-2014 Pieter Kempeneers
This file is part of pktools
diff --git a/src/apps/pkreclass.cc b/src/apps/pkreclass.cc
index e800bfc..b51e13c 100644
--- a/src/apps/pkreclass.cc
+++ b/src/apps/pkreclass.cc
@@ -1,6 +1,6 @@
/**********************************************************************
pkreclass.cc: program to replace pixel values in raster image
-Copyright (C) 2008-2012 Pieter Kempeneers
+Copyright (C) 2008-2014 Pieter Kempeneers
This file is part of pktools
diff --git a/src/apps/pkregann.cc b/src/apps/pkregann.cc
index 1e052a1..cd4ebae 100644
--- a/src/apps/pkregann.cc
+++ b/src/apps/pkregann.cc
@@ -1,6 +1,6 @@
/**********************************************************************
-pkregression_nn.cc: regression with artificial neural network (multi-layer perceptron)
-Copyright (C) 2008-2013 Pieter Kempeneers
+pkregann.cc: regression with artificial neural network (multi-layer perceptron)
+Copyright (C) 2008-2014 Pieter Kempeneers
This file is part of pktools
diff --git a/src/apps/pksetmask.cc b/src/apps/pksetmask.cc
index 6d21e60..cadf101 100644
--- a/src/apps/pksetmask.cc
+++ b/src/apps/pksetmask.cc
@@ -1,6 +1,6 @@
/**********************************************************************
pksetmask.cc: program to apply mask image (set invalid values) to raster image
-Copyright (C) 2008-2012 Pieter Kempeneers
+Copyright (C) 2008-2014 Pieter Kempeneers
This file is part of pktools
diff --git a/src/apps/pksieve.cc b/src/apps/pksieve.cc
index 7d373a5..13673a3 100644
--- a/src/apps/pksieve.cc
+++ b/src/apps/pksieve.cc
@@ -1,6 +1,6 @@
/**********************************************************************
pksieve.cc: program to sieve filter raster image
-Copyright (C) 2008-2012 Pieter Kempeneers
+Copyright (C) 2008-2014 Pieter Kempeneers
This file is part of pktools
diff --git a/src/apps/pkstatascii.cc b/src/apps/pkstatascii.cc
index 6dbb3f1..6dd1d25 100644
--- a/src/apps/pkstatascii.cc
+++ b/src/apps/pkstatascii.cc
@@ -1,6 +1,6 @@
/**********************************************************************
-pkstat.cc: program to calculate basic statistics from raster image
-Copyright (C) 2008-2012 Pieter Kempeneers
+pkstatascii.cc: program to calculate basic statistics from raster image
+Copyright (C) 2008-2014 Pieter Kempeneers
This file is part of pktools
diff --git a/src/apps/pksvm.cc b/src/apps/pksvm.cc
index b428f95..148c0d6 100644
--- a/src/apps/pksvm.cc
+++ b/src/apps/pksvm.cc
@@ -1,6 +1,6 @@
/**********************************************************************
-pkclassify_svm.cc: classify raster image using Support Vector Machine
-Copyright (C) 2008-2012 Pieter Kempeneers
+pksvm.cc: classify raster image using Support Vector Machine
+Copyright (C) 2008-2014 Pieter Kempeneers
This file is part of pktools
@@ -49,13 +49,14 @@ int main(int argc, char *argv[])
//--------------------------- command line options ------------------------------------
Optionpk<string> input_opt("i", "input", "input image");
- Optionpk<string> training_opt("t", "training", "training shape file. A single shape file contains all training features (must be set as: B0, B1, B2,...) for all classes (class numbers identified by label option). Use multiple training files for bootstrap aggregation (alternative to the bag and bsize options, where a random subset is taken from a single training file)");
- Optionpk<string> label_opt("label", "label", "identifier for class label in training shape file.","label");
+ Optionpk<string> training_opt("t", "training", "training vector file. A single vector file contains all training features (must be set as: B0, B1, B2,...) for all classes (class numbers identified by label option). Use multiple training files for bootstrap aggregation (alternative to the bag and bsize options, where a random subset is taken from a single training file)");
+ Optionpk<string> tlayer_opt("tln", "tln", "training layer name(s)");
+ Optionpk<string> label_opt("label", "label", "identifier for class label in training vector file.","label");
Optionpk<unsigned int> balance_opt("bal", "balance", "balance the input data to this number of samples for each class", 0);
Optionpk<bool> random_opt("random", "random", "in case of balance, randomize input data", true);
Optionpk<int> minSize_opt("min", "min", "if number of training pixels is less then min, do not take this class into account (0: consider all classes)", 0);
- Optionpk<double> start_opt("s", "start", "start band sequence number (set to 0)",0);
- Optionpk<double> end_opt("e", "end", "end band sequence number (set to 0 for all bands)", 0);
+ Optionpk<double> start_opt("s", "start", "start band sequence number",0);
+ Optionpk<double> end_opt("e", "end", "end band sequence number (set to 0 to include all bands)", 0);
Optionpk<short> band_opt("b", "band", "band index (starting from 0, either use band option or use start to end)");
Optionpk<double> offset_opt("\0", "offset", "offset value for each spectral band input features: refl[band]=(DN[band]-offset[band])/scale[band]", 0.0);
Optionpk<double> scale_opt("\0", "scale", "scale value for each spectral band input features: refl=(DN[band]-offset[band])/scale[band] (use 0 if scale min and max in each band to -1.0 and 1.0)", 0.0);
@@ -65,9 +66,9 @@ int main(int argc, char *argv[])
Optionpk<std::string> svm_type_opt("svmt", "svmtype", "type of SVM (C_SVC, nu_SVC,one_class, epsilon_SVR, nu_SVR)","C_SVC");
Optionpk<std::string> kernel_type_opt("kt", "kerneltype", "type of kernel function (linear,polynomial,radial,sigmoid) ","radial");
Optionpk<unsigned short> kernel_degree_opt("kd", "kd", "degree in kernel function",3);
- Optionpk<float> gamma_opt("g", "gamma", "gamma in kernel function",0);
+ Optionpk<float> gamma_opt("g", "gamma", "gamma in kernel function",1.0);
Optionpk<float> coef0_opt("c0", "coef0", "coef0 in kernel function",0);
- Optionpk<float> ccost_opt("cc", "ccost", "the parameter C of C_SVC, epsilon_SVR, and nu_SVR",1);
+ Optionpk<float> ccost_opt("cc", "ccost", "the parameter C of C_SVC, epsilon_SVR, and nu_SVR",1000);
Optionpk<float> nu_opt("nu", "nu", "the parameter nu of nu_SVC, one_class SVM, and nu_SVR",0.5);
Optionpk<float> epsilon_loss_opt("eloss", "eloss", "the epsilon in loss function of epsilon_SVR",0.1);
Optionpk<int> cache_opt("cache", "cache", "cache memory size in MB",100);
@@ -99,6 +100,7 @@ int main(int argc, char *argv[])
try{
doProcess=input_opt.retrieveOption(argc,argv);
training_opt.retrieveOption(argc,argv);
+ tlayer_opt.retrieveOption(argc,argv);
label_opt.retrieveOption(argc,argv);
balance_opt.retrieveOption(argc,argv);
random_opt.retrieveOption(argc,argv);
@@ -174,7 +176,7 @@ int main(int argc, char *argv[])
std::cout << "input filename: " << input_opt[0] << std::endl;
if(mask_opt.size())
std::cout << "mask filename: " << mask_opt[0] << std::endl;
- std::cout << "training shape file: " << std::endl;
+ std::cout << "training vector file: " << std::endl;
for(int ifile=0;ifile<training_opt.size();++ifile)
std::cout << training_opt[ifile] << std::endl;
std::cout << "verbose: " << verbose_opt[0] << std::endl;
@@ -249,13 +251,13 @@ int main(int argc, char *argv[])
trainingMap.clear();
trainingPixels.clear();
if(verbose_opt[0]>=1)
- std::cout << "reading imageShape file " << training_opt[0] << std::endl;
+ std::cout << "reading imageVector file " << training_opt[0] << std::endl;
try{
ImgReaderOgr trainingReaderBag(training_opt[ibag]);
if(band_opt.size())
- totalSamples=trainingReaderBag.readDataImageShape(trainingMap,fields,band_opt,label_opt[0],verbose_opt[0]);
+ totalSamples=trainingReaderBag.readDataImageOgr(trainingMap,fields,band_opt,label_opt[0],tlayer_opt,verbose_opt[0]);
else
- totalSamples=trainingReaderBag.readDataImageShape(trainingMap,fields,start_opt[0],end_opt[0],label_opt[0],verbose_opt[0]);
+ totalSamples=trainingReaderBag.readDataImageOgr(trainingMap,fields,start_opt[0],end_opt[0],label_opt[0],tlayer_opt,verbose_opt[0]);
if(trainingMap.size()<2){
string errorstring="Error: could not read at least two classes from training file, did you provide class labels in training sample (see option label)?";
throw(errorstring);
@@ -402,22 +404,29 @@ int main(int argc, char *argv[])
}
map<string,Vector2d<float> >::iterator mapit=trainingMap.begin();
bool doSort=true;
- while(mapit!=trainingMap.end()){
- nameVector.push_back(mapit->first);
- if(classValueMap.size()){
- //check if name in training is covered by classname_opt (values can not be 0)
- if(classValueMap[mapit->first]>0){
- if(cm.getClassIndex(type2string<short>(classValueMap[mapit->first]))<0)
- cm.pushBackClassName(type2string<short>(classValueMap[mapit->first]),doSort);
- }
- else{
- std::cerr << "Error: names in classname option are not complete, please check names in training vector and make sure classvalue is > 0" << std::endl;
- exit(1);
+ try{
+ while(mapit!=trainingMap.end()){
+ nameVector.push_back(mapit->first);
+ if(classValueMap.size()){
+ //check if name in training is covered by classname_opt (values can not be 0)
+ if(classValueMap[mapit->first]>0){
+ if(cm.getClassIndex(type2string<short>(classValueMap[mapit->first]))<0){
+ cm.pushBackClassName(type2string<short>(classValueMap[mapit->first]),doSort);
+ }
+ }
+ else{
+ std::cerr << "Error: names in classname option are not complete, please check names in training vector and make sure classvalue is > 0" << std::endl;
+ exit(1);
+ }
}
+ else
+ cm.pushBackClassName(mapit->first,doSort);
+ ++mapit;
}
- else
- cm.pushBackClassName(mapit->first,doSort);
- ++mapit;
+ }
+ catch(BadConversion conversionString){
+ std::cerr << "Error: did you provide class pairs names (-c) and integer values (-r) for each class in training vector?" << std::endl;
+ exit(1);
}
if(classname_opt.empty()){
std::cerr << "Warning: no class name and value pair provided for all " << nclass << " classes, using string2type<int> instead!" << std::endl;
@@ -428,11 +437,11 @@ int main(int argc, char *argv[])
}
}
- if(priors_opt.size()==nameVector.size()){
- std::cerr << "Warning: please check if priors are provided in correct order!!!" << std::endl;
- for(int iclass=0;iclass<nameVector.size();++iclass)
- std::cerr << nameVector[iclass] << " " << priors_opt[iclass] << std::endl;
- }
+ // if(priors_opt.size()==nameVector.size()){
+ // std::cerr << "Warning: please check if priors are provided in correct order!!!" << std::endl;
+ // for(int iclass=0;iclass<nameVector.size();++iclass)
+ // std::cerr << nameVector[iclass] << " " << priors_opt[iclass] << std::endl;
+ // }
}//if(!ibag)
//Calculate features of training set
@@ -564,16 +573,16 @@ int main(int argc, char *argv[])
if(!verbose_opt[0])
pfnProgress(progress,pszMessage,pProgressArg);
//-------------------------------- open image file ------------------------------------
- bool refIsRaster=false;
+ bool inputIsRaster=false;
ImgReaderOgr imgReaderOgr;
try{
imgReaderOgr.open(input_opt[0]);
imgReaderOgr.close();
}
catch(string errorString){
- refIsRaster=true;
+ inputIsRaster=true;
}
- if(refIsRaster){
+ if(inputIsRaster){
ImgReaderGdal testImage;
try{
if(verbose_opt[0]>=1)
@@ -960,9 +969,9 @@ int main(int argc, char *argv[])
classImageBag.close();
classImageOut.close();
}
- else{//classify shape file
+ else{//classify vector file
cm.clearResults();
- //notice that fields have already been set by readDataImageShape (taking into account appropriate bands)
+ //notice that fields have already been set by readDataImageOgr (taking into account appropriate bands)
for(int ivalidation=0;ivalidation<input_opt.size();++ivalidation){
if(output_opt.size())
assert(output_opt.size()==input_opt.size());
@@ -975,18 +984,20 @@ int main(int argc, char *argv[])
if(verbose_opt[0])
std::cout << "opening img writer and copying fields from img reader" << output_opt[ivalidation] << std::endl;
imgWriterOgr.open(output_opt[ivalidation],imgReaderOgr);
- if(verbose_opt[0])
- std::cout << "creating field class" << std::endl;
- if(classValueMap.size())
- imgWriterOgr.createField("class",OFTInteger);
- else
- imgWriterOgr.createField("class",OFTString);
}
if(verbose_opt[0])
cout << "number of layers in input ogr file: " << imgReaderOgr.getLayerCount() << endl;
for(int ilayer=0;ilayer<imgReaderOgr.getLayerCount();++ilayer){
if(verbose_opt[0])
cout << "processing input layer " << ilayer << endl;
+ if(output_opt.size()){
+ if(verbose_opt[0])
+ std::cout << "creating field class" << std::endl;
+ if(classValueMap.size())
+ imgWriterOgr.createField("class",OFTInteger,ilayer);
+ else
+ imgWriterOgr.createField("class",OFTString,ilayer);
+ }
unsigned int nFeatures=imgReaderOgr.getFeatureCount(ilayer);
unsigned int ifeature=0;
progress=0;
@@ -1001,11 +1012,11 @@ int main(int argc, char *argv[])
}
OGRFeature *poDstFeature = NULL;
if(output_opt.size()){
- poDstFeature=imgWriterOgr.createFeature();
+ poDstFeature=imgWriterOgr.createFeature(ilayer);
if( poDstFeature->SetFrom( poFeature, TRUE ) != OGRERR_NONE ){
CPLError( CE_Failure, CPLE_AppDefined,
"Unable to translate feature %d from layer %s.\n",
- poFeature->GetFID(), imgWriterOgr.getLayerName().c_str() );
+ poFeature->GetFID(), imgWriterOgr.getLayerName(ilayer).c_str() );
OGRFeature::DestroyFeature( poFeature );
OGRFeature::DestroyFeature( poDstFeature );
}
@@ -1113,10 +1124,10 @@ int main(int argc, char *argv[])
}
CPLErrorReset();
if(output_opt.size()){
- if(imgWriterOgr.createFeature( poDstFeature ) != OGRERR_NONE){
+ if(imgWriterOgr.createFeature(poDstFeature,ilayer) != OGRERR_NONE){
CPLError( CE_Failure, CPLE_AppDefined,
"Unable to translate feature %d from layer %s.\n",
- poFeature->GetFID(), imgWriterOgr.getLayerName().c_str() );
+ poFeature->GetFID(), imgWriterOgr.getLayerName(ilayer).c_str() );
OGRFeature::DestroyFeature( poDstFeature );
OGRFeature::DestroyFeature( poDstFeature );
}
diff --git a/src/base/Makefile.am b/src/base/Makefile.am
index d641e68..d2228de 100644
--- a/src/base/Makefile.am
+++ b/src/base/Makefile.am
@@ -27,7 +27,7 @@ libbase_ladir = $(includedir)/pktools/base
libbase_la_LDFLAGS = -version-info $(PKTOOLS_SO_VERSION) $(AM_LDFLAGS)
# the list of header files that belong to the library (to be installed later)
-libbase_la_HEADERS = $(top_srcdir)/config.h Vector2d.h IndexValue.h Optionpk.h PosValue.h
+libbase_la_HEADERS = Vector2d.h IndexValue.h Optionpk.h PosValue.h
# the sources to add to the library and to add to the source distribution
###############################################################################
diff --git a/src/base/Makefile.in b/src/base/Makefile.in
index f985404..c716975 100644
--- a/src/base/Makefile.in
+++ b/src/base/Makefile.in
@@ -301,7 +301,7 @@ libbase_ladir = $(includedir)/pktools/base
libbase_la_LDFLAGS = -version-info $(PKTOOLS_SO_VERSION) $(AM_LDFLAGS)
# the list of header files that belong to the library (to be installed later)
-libbase_la_HEADERS = $(top_srcdir)/config.h Vector2d.h IndexValue.h Optionpk.h PosValue.h
+libbase_la_HEADERS = Vector2d.h IndexValue.h Optionpk.h PosValue.h
# the sources to add to the library and to add to the source distribution
###############################################################################
diff --git a/src/base/Vector2d.h b/src/base/Vector2d.h
index 12ac813..bbbcafa 100644
--- a/src/base/Vector2d.h
+++ b/src/base/Vector2d.h
@@ -49,6 +49,7 @@ public:
void selectCol(int col, T* output) const;
std::vector<T> selectCol(int col);
void selectCols(const std::list<int> &cols, Vector2d<T> &output) const;
+ void setMask(const Vector2d<T> &mask, T msknodata, T nodata=0);
void transpose(Vector2d<T> &output) const{
output.resize(nCols(),nRows());
for(int irow=0;irow<nRows();++irow){
@@ -62,6 +63,7 @@ public:
void scale(const std::vector<double> &scaleVector, const std::vector<double> &offsetVector, Vector2d<T>& scaledOutput);
void scale(const T lbound, const T ubound, std::vector<double> &scaleVector, std::vector<double> &offsetVector, Vector2d<T>& scaledOutput);
Vector2d<T> operator=(const Vector2d<T>& v1);
+ Vector2d<T> operator+=(const Vector2d<T>& v1);
// std::ostream& operator<<(std::ostream& os, const Vector2d<T>& v);
// template<class T> std::ostream& operator<<(std::ostream& os, const Vector2d<T>& v);
template<class T1> friend std::ostream& operator<<(std::ostream & os, const Vector2d<T1>& v);
@@ -99,6 +101,15 @@ template<class T> Vector2d<T> Vector2d<T>::operator=(const Vector2d<T>& v1){
}
}
+template<class T> Vector2d<T> Vector2d<T>::operator+=(const Vector2d<T>& v1){
+ assert(v1.nRows()==nRows());
+ assert(v1.nCols()==nCols());
+ for(int irow=0;irow<nRows();++irow)
+ for(int icol=0;icol<nCols();++icol)
+ (*this)[irow][icol]+=v1[irow][icol];
+ return *this;
+}
+
template<class T> Vector2d<T>::Vector2d(int nrow)
: std::vector< std::vector<T> >(nrow)
{
@@ -192,6 +203,16 @@ template<class T> void Vector2d<T>::selectCols(const std::list<int> &cols)
(*this)[irow].erase(((*this)[irow]).begin()+icol);
}
+template<class T> void Vector2d<T>::setMask(const Vector2d<T> &mask, T msknodata, T nodata)
+{
+ assert(mask.nRows()==nRows());
+ assert(mask.nCols()==nCols());
+ for(int irow=0;irow<this->size();++irow)
+ for(int icol=0;icol<((*this)[irow]).size()-1;++icol)
+ if(mask[irow][icol]==msknodata)
+ (*this)[irow][icol]=nodata;
+}
+
template<class T1> std::ostream& operator<<(std::ostream& os, const Vector2d<T1>& v)
{
for(int irow=0;irow<v.size();++irow){
diff --git a/src/imageclasses/ImgReaderGdal.cc b/src/imageclasses/ImgReaderGdal.cc
index 6bd8093..63d25bd 100644
--- a/src/imageclasses/ImgReaderGdal.cc
+++ b/src/imageclasses/ImgReaderGdal.cc
@@ -258,7 +258,7 @@ bool ImgReaderGdal::getBoundingBox(double& ulx, double& uly, double& lrx, double
}
}
-bool ImgReaderGdal::getCentrePos(double& x, double& y) const
+bool ImgReaderGdal::getCenterPos(double& x, double& y) const
{
double gt[6];// { 444720, 30, 0, 3751320, 0, -30 };
m_gds->GetGeoTransform(gt);
@@ -318,7 +318,7 @@ bool ImgReaderGdal::geo2image(double x, double y, double& i, double& j) const
}
}
-//x and y represent centre of pixel, return true if image is georeferenced
+//x and y represent center of pixel, return true if image is georeferenced
bool ImgReaderGdal::image2geo(double i, double j, double& x, double& y) const
{
double gt[6];// { 444720, 30, 0, 3751320, 0, -30 };
@@ -619,15 +619,15 @@ void ImgReaderGdal::getRefPix(double& refX, double &refY, int band) const
}
}
if(isGeoRef()){
- //reference coordinate is lower left corner of pixel in centre of gravity
+ //reference coordinate is lower left corner of pixel in center of gravity
//we need geo coordinates for exactly this location: validCol(Row)/nvalidCol(Row)-0.5
double cgravi=validCol/nvalidCol-0.5;
double cgravj=validRow/nvalidRow-0.5;
double refpixeli=floor(cgravi);
double refpixelj=ceil(cgravj-1);
- //but image2geo provides location at centre of pixel (shifted half pixel right down)
+ //but image2geo provides location at center of pixel (shifted half pixel right down)
image2geo(refpixeli,refpixelj,refX,refY);
- //refX and refY now refer to centre of gravity pixel
+ //refX and refY now refer to center of gravity pixel
refX-=0.5*getDeltaX();//shift to left corner
refY-=0.5*getDeltaY();//shift to lower left corner
}
diff --git a/src/imageclasses/ImgReaderGdal.h b/src/imageclasses/ImgReaderGdal.h
index c3303c1..2d941a6 100644
--- a/src/imageclasses/ImgReaderGdal.h
+++ b/src/imageclasses/ImgReaderGdal.h
@@ -51,7 +51,7 @@ public:
std::string getMetadataItem() const;
std::string getImageDescription() const;
bool getBoundingBox (double& ulx, double& uly, double& lrx, double& lry) const;
- bool getCentrePos(double& x, double& y) const;
+ bool getCenterPos(double& x, double& y) const;
double getUlx() const {double ulx, uly, lrx,lry;getBoundingBox(ulx,uly,lrx,lry);return(ulx);};
double getUly() const {double ulx, uly, lrx,lry;getBoundingBox(ulx,uly,lrx,lry);return(uly);};
double getLrx() const {double ulx, uly, lrx,lry;getBoundingBox(ulx,uly,lrx,lry);return(lrx);};
diff --git a/src/imageclasses/ImgReaderOgr.cc b/src/imageclasses/ImgReaderOgr.cc
index 500dac9..db5d7cc 100644
--- a/src/imageclasses/ImgReaderOgr.cc
+++ b/src/imageclasses/ImgReaderOgr.cc
@@ -217,11 +217,12 @@ std::ostream& operator<<(std::ostream& theOstream, ImgReaderOgr& theImageReader)
// }
// }
-unsigned int ImgReaderOgr::readDataImageShape(std::map<std::string,Vector2d<float> > &mapPixels, //[classNr][pixelNr][bandNr],
- std::vector<std::string>& fields,
- const std::vector<short>& bands,
- const std::string& label,
- int verbose)
+unsigned int ImgReaderOgr::readDataImageOgr(std::map<std::string,Vector2d<float> > &mapPixels, //[classNr][pixelNr][bandNr],
+ std::vector<std::string>& fields,
+ const std::vector<short>& bands,
+ const std::string& label,
+ const std::vector<std::string>& layers,
+ int verbose)
{
mapPixels.clear();
int nsample=0;
@@ -230,6 +231,10 @@ unsigned int ImgReaderOgr::readDataImageShape(std::map<std::string,Vector2d<floa
if(verbose)
std::cout << "reading shape file " << m_filename << std::endl;
for(int ilayer=0;ilayer<getLayerCount();++ilayer){
+ std::string currentLayername=getLayer(ilayer)->GetName();
+ if(layers.size())
+ if(find(layers.begin(),layers.end(),currentLayername)==layers.end())
+ continue;
try{
//only retain bands in fields
getFields(fields,ilayer);
@@ -295,12 +300,13 @@ unsigned int ImgReaderOgr::readDataImageShape(std::map<std::string,Vector2d<floa
return totalSamples;
}
-unsigned int ImgReaderOgr::readDataImageShape(std::map<std::string,Vector2d<float> > &mapPixels, //[classNr][pixelNr][bandNr],
- std::vector<std::string>& fields,
- double start,
- double end,
- const std::string& label,
- int verbose)
+unsigned int ImgReaderOgr::readDataImageOgr(std::map<std::string,Vector2d<float> > &mapPixels, //[classNr][pixelNr][bandNr],
+ std::vector<std::string>& fields,
+ double start,
+ double end,
+ const std::string& label,
+ const std::vector<std::string>& layers,
+ int verbose)
{
mapPixels.clear();
int nsample=0;
@@ -309,6 +315,10 @@ unsigned int ImgReaderOgr::readDataImageShape(std::map<std::string,Vector2d<floa
if(verbose)
std::cout << "reading shape file " << m_filename << std::endl;
for(int ilayer=0;ilayer<getLayerCount();++ilayer){
+ std::string currentLayername=getLayer(ilayer)->GetName();
+ if(layers.size())
+ if(find(layers.begin(),layers.end(),currentLayername)==layers.end())
+ continue;
try{
//only retain bands in fields
getFields(fields,ilayer);
diff --git a/src/imageclasses/ImgReaderOgr.h b/src/imageclasses/ImgReaderOgr.h
index 36653e6..f15f905 100644
--- a/src/imageclasses/ImgReaderOgr.h
+++ b/src/imageclasses/ImgReaderOgr.h
@@ -48,18 +48,20 @@ public:
template <typename T> int readData(Vector2d<T>& data, const OGRFieldType& fieldType, std::vector<std::string>& fields, int layer=0, bool pos=false, bool verbose=false);//default layer 0 and no pos information in data
template <typename T> int readData(std::map<int,Vector2d<T> >& data, const OGRFieldType& fieldType, std::vector<std::string>& fields, const std::string& label, int layer=0, bool pos=false, bool verbose=false);//default layer 0 and no pos information in data
template <typename T> int readData(std::map<std::string,Vector2d<T> >& data, const OGRFieldType& fieldType, std::vector<std::string>& fields, const std::string& label, int layer=0, bool pos=false, bool verbose=false);//default layer 0 and no pos information in data
- unsigned int readDataImageShape(std::map<std::string,Vector2d<float> > &mapPixels, //[classNr][pixelNr][bandNr],
- std::vector<std::string>& fields,
- const std::vector<short>& bands,
- const std::string& label,
- int verbose=false);
+ unsigned int readDataImageOgr(std::map<std::string,Vector2d<float> > &mapPixels, //[classNr][pixelNr][bandNr],
+ std::vector<std::string>& fields,
+ const std::vector<short>& bands,
+ const std::string& label,
+ const std::vector<std::string>& layers,
+ int verbose=false);
- unsigned int readDataImageShape(std::map<std::string,Vector2d<float> > &mapPixels, //[classNr][pixelNr][bandNr],
- std::vector<std::string>& fields,
- double start,
- double end,
- const std::string& label,
- int verbose=false);
+ unsigned int readDataImageOgr(std::map<std::string,Vector2d<float> > &mapPixels, //[classNr][pixelNr][bandNr],
+ std::vector<std::string>& fields,
+ double start,
+ double end,
+ const std::string& label,
+ const std::vector<std::string>& layers,
+ int verbose=false);
void shape2ascii(std::ostream& theOstream, const std::string& pointname, int layer=0, bool verbose=false);
unsigned long int getFeatureCount(int layer=0) const;
diff --git a/src/imageclasses/ImgWriterOgr.cc b/src/imageclasses/ImgWriterOgr.cc
index d5e5de2..4f040cf 100644
--- a/src/imageclasses/ImgWriterOgr.cc
+++ b/src/imageclasses/ImgWriterOgr.cc
@@ -45,7 +45,7 @@ ImgWriterOgr::ImgWriterOgr(const std::string& filename, ImgReaderOgr& imgReaderO
for(int ilayer=0;ilayer<nlayer;++ilayer){
std::string layername = imgReaderOgr.getLayer(ilayer)->GetName();
createLayer(layername,imgReaderOgr.getProjection(),imgReaderOgr.getGeometryType(),NULL);
- copyFields(imgReaderOgr,ilayer);
+ copyFields(imgReaderOgr,ilayer,ilayer);
}
}
@@ -58,7 +58,7 @@ ImgWriterOgr::ImgWriterOgr(const std::string& filename, ImgReaderOgr& imgReaderO
for(int ilayer=0;ilayer<nlayer;++ilayer){
std::string layername = imgReaderOgr.getLayer(ilayer)->GetName();
createLayer(layername,imgReaderOgr.getProjection(),imgReaderOgr.getGeometryType(),NULL);
- copyFields(imgReaderOgr,ilayer);
+ copyFields(imgReaderOgr,ilayer,ilayer);
if(copyData){
OGRFeature *poFeature;
while(true){// (poFeature = imgReaderOgr.getLayer()->GetNextFeature()) != NULL ){
@@ -112,7 +112,7 @@ void ImgWriterOgr::open(const std::string& filename, ImgReaderOgr& imgReaderOgr)
for(int ilayer=0;ilayer<nlayer;++ilayer){
std::string layername = imgReaderOgr.getLayer(ilayer)->GetName();
createLayer(layername,imgReaderOgr.getProjection(),imgReaderOgr.getGeometryType(),NULL);
- copyFields(imgReaderOgr,ilayer);
+ copyFields(imgReaderOgr,ilayer,ilayer);
}
}
@@ -262,16 +262,15 @@ int ImgWriterOgr::getFields(std::vector<OGRFieldDefn*>& fields, int layer) const
return(fields.size());
}
-void ImgWriterOgr::copyFields(const ImgReaderOgr& imgReaderOgr, int theLayer){
- if(theLayer<0)
- theLayer=m_datasource->GetLayerCount()-1;//get back layer
+void ImgWriterOgr::copyFields(const ImgReaderOgr& imgReaderOgr, int srcLayer, int targetLayer){
+ if(targetLayer<0)
+ targetLayer=m_datasource->GetLayerCount()-1;//get back layer
//get fields from imgReaderOgr
std::vector<OGRFieldDefn*> fields;
- imgReaderOgr.getFields(fields,theLayer);
-// OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
+ imgReaderOgr.getFields(fields,srcLayer);
for(int iField=0;iField<fields.size();++iField){
- if(m_datasource->GetLayer(theLayer)->CreateField(fields[iField]) != OGRERR_NONE ){
+ if(m_datasource->GetLayer(targetLayer)->CreateField(fields[iField]) != OGRERR_NONE ){
std::ostringstream es;
es << "Creating field " << fields[iField]->GetNameRef() << " failed";
std::string errorString=es.str();
@@ -602,7 +601,7 @@ int ImgWriterOgr::addData(const ImgReaderGdal& imgReader, int theLayer, bool ver
for(int iband=0;iband<imgReader.nrOfBand();++iband){
std::ostringstream fs;
fs << "band" << iband;
- createField(fs.str(),OFTReal);
+ createField(fs.str(),OFTReal,theLayer);
}
OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
int nfield=poFDefn->GetFieldCount();
diff --git a/src/imageclasses/ImgWriterOgr.h b/src/imageclasses/ImgWriterOgr.h
index 6d54a69..327b9a7 100644
--- a/src/imageclasses/ImgWriterOgr.h
+++ b/src/imageclasses/ImgWriterOgr.h
@@ -51,7 +51,7 @@ public:
std::string getLayerName(int layer=0){return m_datasource->GetLayer(layer)->GetLayerDefn()->GetName();};
int getFields(std::vector<std::string>& fields, int layer=0) const;
int getFields(std::vector<OGRFieldDefn*>& fields, int layer=0) const;
- void copyFields(const ImgReaderOgr& imgReaderOgr, int theLayer=0);//default: get back layer
+ void copyFields(const ImgReaderOgr& imgReaderOgr, int srcLayer=0, int targetLayer=0);//default: get back layer
void addLineString(std::vector<OGRPoint*>& points, const std::string& fieldName, const std::string& theId, int layer=0);
void addRing(std::vector<OGRPoint*>& points, const std::string& fieldName, int theId, int layer=0);
void addLineString(std::vector<OGRPoint*>& points, const std::string& fieldName, int theId, int layer=0);
diff --git a/src/lasclasses/FileReaderLas.cc b/src/lasclasses/FileReaderLas.cc
index aceaac1..f6873c0 100644
--- a/src/lasclasses/FileReaderLas.cc
+++ b/src/lasclasses/FileReaderLas.cc
@@ -83,9 +83,11 @@ void FileReaderLas::close(void)
void FileReaderLas::setCodec(const std::string& filename){
m_ifstream = new(std::ifstream);
m_ifstream->open(m_filename.c_str(),std::ios::in|std::ios::binary);
- m_reader = new liblas::Reader(*m_ifstream);
+ liblas::ReaderFactory f;
+ liblas::Reader reader = f.CreateWithStream(*m_ifstream);
+ m_reader=new liblas::Reader(reader);
+ // m_reader = new liblas::Reader(*m_ifstream);
//Note: It is possible to use the basic liblas::Reader constructor that takes in a std::istream, but it will not be able to account for the fact that the file might be compressed. Using the ReaderFactory will take care of all of this for you.
- // liblas::ReaderFactory rfactory;
// m_reader=&rfactory.CreateWithStream(ifs);
}
--
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