[dans-gdal-scripts] 01/05: Imported Upstream version 0.24
Bas Couwenberg
sebastic at debian.org
Sat Oct 15 18:56:52 UTC 2016
This is an automated email from the git hooks/post-receive script.
sebastic pushed a commit to branch master
in repository dans-gdal-scripts.
commit 260ca371a8110b21e958751a58010b255df627ba
Author: Bas Couwenberg <sebastic at xs4all.nl>
Date: Sat Oct 15 20:39:34 2016 +0200
Imported Upstream version 0.24
---
ChangeLog | 8 +
TODO | 24 +-
configure.ac | 3 +-
m4/boost.m4 | 482 ++++++++++++++++-----
src/Makefile.am | 12 +-
src/common.cc | 4 +
src/common.h | 4 +
src/datatype_conversion.cc | 65 +++
src/datatype_conversion.h | 49 +++
src/default_palette.h | 33 ++
src/dp.cc | 155 ++++---
src/dp.h | 7 +
src/gdal_contrast_stretch.cc | 32 +-
src/gdal_dem2rgb.cc | 16 +-
src/gdal_get_projected_bounds.cc | 2 +-
src/gdal_landsat_pansharp.cc | 2 +
src/gdal_merge_simple.cc | 2 +
src/gdal_merge_vrt.cc | 2 +
src/gdal_raw2geotiff.cc | 2 +
src/gdal_trace_outline.cc | 131 +++---
src/mask-tracer.cc | 68 +--
src/mask.cc | 273 ++++--------
src/mask.h | 72 +--
src/ndv.cc | 160 +++----
src/ndv.h | 65 ++-
src/polygon.cc | 63 ++-
src/polygon.h | 160 +++++++
src/raster_features.cc | 440 +++++++++++++++++++
src/raster_features.h | 118 +++++
src/tests/.gitignore | 4 +
src/tests/good_test1_1.ppm | Bin 0 -> 198162 bytes
src/tests/{good_test2_2.ppm => good_test1_2.ppm} | Bin 198162 -> 198162 bytes
src/tests/good_test1_3.ppm | Bin 0 -> 198162 bytes
src/tests/good_test1_3_classify.dbf | Bin 0 -> 284 bytes
src/tests/good_test1_3_classify.shp | Bin 0 -> 75220 bytes
src/tests/good_test1_3_classify.shx | Bin 0 -> 140 bytes
src/tests/good_test1_3_classify_pal.dbf | Bin 0 -> 448 bytes
src/tests/good_test1_3_classify_pal.shp | Bin 0 -> 75220 bytes
src/tests/good_test1_3_classify_pal.shx | Bin 0 -> 140 bytes
src/tests/good_test1_4-rect.ppm | Bin 1277157 -> 1277157 bytes
src/tests/{good_test2_4.ppm => good_test1_4.ppm} | Bin 1277157 -> 1277157 bytes
src/tests/{good_test2_5.ppm => good_test1_5.ppm} | Bin 1177863 -> 1177863 bytes
src/tests/good_test1_double.dbf | Bin 0 -> 390 bytes
src/tests/good_test1_double.shp | Bin 0 -> 464856 bytes
src/tests/good_test1_double.shx | Bin 0 -> 204 bytes
src/tests/good_test1_double.wkt | 13 +
src/tests/good_test1_double_clip.wkt | 3 +
src/tests/good_test1_features.dbf | Bin 0 -> 749 bytes
src/tests/good_test1_features.shp | Bin 0 -> 2009632 bytes
src/tests/good_test1_features.shx | Bin 0 -> 260 bytes
src/tests/good_test1_features.wkt | 20 +
src/tests/good_test1_gradient_ndv.pbm | Bin 0 -> 8203 bytes
src/tests/good_test1_gradient_ndv_inv.pbm | Bin 0 -> 8203 bytes
src/tests/good_test1_gradient_valid.pbm | Bin 0 -> 8203 bytes
src/tests/good_test1_maze.ppm | Bin 0 -> 771138 bytes
src/tests/good_test1_nan1.pbm | Bin 0 -> 8203 bytes
src/tests/good_test1_nan2.pbm | Bin 0 -> 8203 bytes
src/tests/good_test1_nan3.pbm | Bin 0 -> 8203 bytes
src/tests/good_test1_nan4.pbm | Bin 0 -> 8203 bytes
src/tests/good_test1_noise.ppm | Bin 0 -> 771138 bytes
src/tests/good_test1_noise_dp3.ppm | Bin 0 -> 771138 bytes
src/tests/good_test2_1.ppm | Bin 198162 -> 198162 bytes
src/tests/good_test2_2.ppm | Bin 198162 -> 198162 bytes
src/tests/good_test2_3.ppm | Bin 198162 -> 198162 bytes
src/tests/good_test2_4.ppm | Bin 1277157 -> 1277157 bytes
src/tests/good_test2_5.ppm | Bin 1177863 -> 1177863 bytes
src/tests/good_test2_maze.ppm | Bin 771138 -> 771138 bytes
src/tests/good_test2_noise.ppm | Bin 771138 -> 771138 bytes
src/tests/good_test2_noise_dp3.ppm | Bin 771138 -> 771138 bytes
src/tests/good_test3_dem.tif | Bin 0 -> 618214 bytes
.../{good_test2_4.ppm => good_test3_histeq.tif} | Bin 1277157 -> 1699534 bytes
src/tests/test1.sh | 94 +++-
src/tests/test2.sh | 8 +-
src/tests/test3.sh | 8 +-
src/tests/{pal.tif => testcase_3_paletted.tif} | Bin
src/tests/testcase_double.tif | Bin 0 -> 824313 bytes
src/tests/testcase_features.png | Bin 0 -> 85946 bytes
77 files changed, 1936 insertions(+), 668 deletions(-)
diff --git a/ChangeLog b/ChangeLog
index 62ec57d..5fd6fab 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,11 @@
+=== Version 0.24
+ gdal_trace_outline no longer hangs at the "fixing topology" stage
+ gdal_trace_outline -classify can now handle arbitrary datatypes, multiple bands, and ndv
+ bugfix: "-ndv '1 1 1' -ndv '2 2 2'" would also match '1 1 2' values, for instance
+ ndv option accepts '*' as alias for '-Inf..Inf'
+ fixed boost version detection
+ various minor improvements
+
=== Version 0.23
Fix for compiler warnings/errors.
Fix finding of GDAL headers.
diff --git a/TODO b/TODO
index b04531c..98fa973 100644
--- a/TODO
+++ b/TODO
@@ -1,13 +1,19 @@
-gdal_list_corners fails on venus gtdr.tif (lat>90)
-dem2rgb and others: maybe NDV should be independent of band selection
+gdal_contrast_stretch:
+ * allow 16-bit output
+ * several inputs for contrast stretch
+ * manual min/max for contrast stretch
+ (gdal_translate -scale already does this)
+
+gdal_trace_outline:
+ * for feature classification, trace all feature classes simultaneously (for
+ much better speed)
+ * run erosion several times for outline tracer
+ * expose options for fuzzy rectangle bounds finder
+ * use concave hull instead of the current excursions pincher
+ * outline tracer should call OGR_G_IsValid on result
+
+gdal_list_corners: fails on venus gtdr.tif (lat>90)
dem2rgb: prominant warning if resolution not known
dem2rgb: should be operable with resolution in absence of origin
-run erosion several times for outline tracer
gdal_merge_simple: usage text wrong?
-beveller is inefficient (maybe O(N^2))
don't use %f, use %g with the right number of digits
-several inputs for contrast stretch
-manual min/max for contrast stretch
-better error message for NDV>255 in contrast stretch
-expose options for fuzzy rectangle bounds finder
-outline tracer should call OGR_G_IsValid on result
diff --git a/configure.ac b/configure.ac
index df97915..039ee26 100644
--- a/configure.ac
+++ b/configure.ac
@@ -24,7 +24,7 @@
AC_PREREQ(2.37)
-AC_INIT([dans-gdal-scripts], [0.23], [dan at stahlke.org])
+AC_INIT([dans-gdal-scripts], [0.24], [dan at stahlke.org])
AC_CONFIG_SRCDIR([src/common.h])
AC_CONFIG_HEADER([config.h])
@@ -89,7 +89,6 @@ AC_FUNC_STRTOD
AC_CHECK_FUNCS([memset sqrt strtol])
AC_CHECK_LIB([m], [pow], , AC_MSG_ERROR([math library is required]))
-AC_CHECK_LIB([proj], [pj_init], , AC_MSG_ERROR([proj4 library is required]))
AC_CONFIG_FILES([Makefile
src/Makefile])
diff --git a/m4/boost.m4 b/m4/boost.m4
index 3d4e47c..bd0a9bd 100644
--- a/m4/boost.m4
+++ b/m4/boost.m4
@@ -1,5 +1,5 @@
# boost.m4: Locate Boost headers and libraries for autoconf-based projects.
-# Copyright (C) 2007, 2008, 2009, 2010, 2011 Benoit Sigoure <tsuna at lrde.epita.fr>
+# Copyright (C) 2007-2011, 2014 Benoit Sigoure <tsuna at lrde.epita.fr>
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
@@ -22,7 +22,7 @@
# along with this program. If not, see <http://www.gnu.org/licenses/>.
m4_define([_BOOST_SERIAL], [m4_translit([
-# serial 16
+# serial 24
], [#
], [])])
@@ -59,7 +59,8 @@ m4_pattern_forbid([^_?(BOOST|Boost)_])
# It could be useful to turn this into a macro which extracts the
# value of any macro.
m4_define([_BOOST_SED_CPP],
-[AC_LANG_PREPROC_REQUIRE()dnl
+[AC_LANG_PUSH([C++])dnl
+AC_LANG_PREPROC_REQUIRE()dnl
AC_REQUIRE([AC_PROG_SED])dnl
AC_LANG_CONFTEST([AC_LANG_SOURCE([[$2]])])
AS_IF([dnl eval is necessary to expand ac_cpp.
@@ -71,13 +72,31 @@ dnl strip `\n' with backquotes, not the `\r'. This results in
dnl boost_cv_lib_version='1_37\r' for instance, which breaks
dnl everything else.
dnl Cannot use 'dnl' after [$4] because a trailing dnl may break AC_CACHE_CHECK
+dnl
+dnl Beware that GCC 5, when expanding macros, may embed # line directives
+dnl a within single line:
+dnl
+dnl # 1 "conftest.cc"
+dnl # 1 "<built-in>"
+dnl # 1 "<command-line>"
+dnl # 1 "conftest.cc"
+dnl # 1 "/opt/local/include/boost/version.hpp" 1 3
+dnl # 2 "conftest.cc" 2
+dnl boost-lib-version =
+dnl # 2 "conftest.cc" 3
+dnl "1_56"
+dnl
+dnl So get rid of the # lines, and glue the remaining ones together.
(eval "$ac_cpp conftest.$ac_ext") 2>&AS_MESSAGE_LOG_FD |
+ grep -v '#' |
tr -d '\r' |
+ tr -s '\n' ' ' |
$SED -n -e "$1" >conftest.i 2>&1],
[$3],
[$4])
rm -rf conftest*
-])# AC_EGREP_CPP
+AC_LANG_POP([C++])dnl
+])# _BOOST_SED_CPP
@@ -206,7 +225,7 @@ AC_LANG_POP([C++])dnl
AC_CACHE_CHECK([for Boost's header version],
[boost_cv_lib_version],
[m4_pattern_allow([^BOOST_LIB_VERSION$])dnl
- _BOOST_SED_CPP([/^boost-lib-version = /{s///;s/\"//g;p;q;}],
+ _BOOST_SED_CPP([[/^boost-lib-version = /{s///;s/[\" ]//g;p;q;}]],
[#include <boost/version.hpp>
boost-lib-version = BOOST_LIB_VERSION],
[boost_cv_lib_version=`cat conftest.i`])])
@@ -214,24 +233,26 @@ boost-lib-version = BOOST_LIB_VERSION],
boost_major_version=`echo "$boost_cv_lib_version" | sed 's/_//;s/_.*//'`
case $boost_major_version in #(
'' | *[[!0-9]]*)
- AC_MSG_ERROR([invalid value: boost_major_version=$boost_major_version])
+ AC_MSG_ERROR([invalid value: boost_major_version='$boost_major_version'])
;;
esac
fi
CPPFLAGS=$boost_save_CPPFLAGS
])# BOOST_REQUIRE
+
# BOOST_STATIC()
# --------------
# Add the "--enable-static-boost" configure argument. If this argument is given
# on the command line, static versions of the libraries will be looked up.
AC_DEFUN([BOOST_STATIC],
[AC_ARG_ENABLE([static-boost],
- [AC_HELP_STRING([--enable-static-boost],
+ [AS_HELP_STRING([--enable-static-boost],
[Prefer the static boost libraries over the shared ones [no]])],
[enable_static_boost=yes],
[enable_static_boost=no])])# BOOST_STATIC
+
# BOOST_FIND_HEADER([HEADER-NAME], [ACTION-IF-NOT-FOUND], [ACTION-IF-FOUND])
# --------------------------------------------------------------------------
# Wrapper around AC_CHECK_HEADER for Boost headers. Useful to check for
@@ -264,14 +285,16 @@ fi
])# BOOST_FIND_HEADER
-# BOOST_FIND_LIB([LIB-NAME], [PREFERRED-RT-OPT], [HEADER-NAME], [CXX-TEST],
-# [CXX-PROLOGUE])
-# -------------------------------------------------------------------------
-# Look for the Boost library LIB-NAME (e.g., LIB-NAME = `thread', for
-# libboost_thread). Check that HEADER-NAME works and check that
-# libboost_LIB-NAME can link with the code CXX-TEST. The optional argument
-# CXX-PROLOGUE can be used to include some C++ code before the `main'
-# function.
+# BOOST_FIND_LIBS([COMPONENT-NAME], [CANDIDATE-LIB-NAMES],
+# [PREFERRED-RT-OPT], [HEADER-NAME], [CXX-TEST],
+# [CXX-PROLOGUE])
+# --------------------------------------------------------------
+# Look for the Boost library COMPONENT-NAME (e.g., `thread', for
+# libboost_thread) under the possible CANDIDATE-LIB-NAMES (e.g.,
+# "thread_win32 thread"). Check that HEADER-NAME works and check that
+# libboost_LIB-NAME can link with the code CXX-TEST. The optional
+# argument CXX-PROLOGUE can be used to include some C++ code before
+# the `main' function.
#
# Invokes BOOST_FIND_HEADER([HEADER-NAME]) (see above).
#
@@ -285,7 +308,7 @@ fi
# builds. Some sample values for PREFERRED-RT-OPT: (nothing), mt, d, mt-d, gdp
# ... If you want to make sure you have a specific version of Boost
# (eg, >= 1.33) you *must* invoke BOOST_REQUIRE before this macro.
-AC_DEFUN([BOOST_FIND_LIB],
+AC_DEFUN([BOOST_FIND_LIBS],
[AC_REQUIRE([BOOST_REQUIRE])dnl
AC_REQUIRE([_BOOST_FIND_COMPILER_TAG])dnl
AC_REQUIRE([BOOST_STATIC])dnl
@@ -299,32 +322,69 @@ AS_VAR_PUSHDEF([Boost_lib], [boost_cv_lib_$1])dnl
AS_VAR_PUSHDEF([Boost_lib_LDFLAGS], [boost_cv_lib_$1_LDFLAGS])dnl
AS_VAR_PUSHDEF([Boost_lib_LDPATH], [boost_cv_lib_$1_LDPATH])dnl
AS_VAR_PUSHDEF([Boost_lib_LIBS], [boost_cv_lib_$1_LIBS])dnl
-BOOST_FIND_HEADER([$3])
+BOOST_FIND_HEADER([$4])
boost_save_CPPFLAGS=$CPPFLAGS
CPPFLAGS="$CPPFLAGS $BOOST_CPPFLAGS"
-# Now let's try to find the library. The algorithm is as follows: first look
-# for a given library name according to the user's PREFERRED-RT-OPT. For each
-# library name, we prefer to use the ones that carry the tag (toolset name).
-# Each library is searched through the various standard paths were Boost is
-# usually installed. If we can't find the standard variants, we try to
-# enforce -mt (for instance on MacOSX, libboost_threads.dylib doesn't exist
-# but there's -obviously- libboost_threads-mt.dylib).
AC_CACHE_CHECK([for the Boost $1 library], [Boost_lib],
- [Boost_lib=no
- case "$2" in #(
- mt | mt-) boost_mt=-mt; boost_rtopt=;; #(
- mt* | mt-*) boost_mt=-mt; boost_rtopt=`expr "X$2" : 'Xmt-*\(.*\)'`;; #(
- *) boost_mt=; boost_rtopt=$2;;
+ [_BOOST_FIND_LIBS($@)])
+case $Boost_lib in #(
+ (no) _AC_MSG_LOG_CONFTEST
+ AC_MSG_ERROR([cannot find the flags to link with Boost $1])
+ ;;
+esac
+AC_SUBST(AS_TR_CPP([BOOST_$1_LDFLAGS]), [$Boost_lib_LDFLAGS])dnl
+AC_SUBST(AS_TR_CPP([BOOST_$1_LDPATH]), [$Boost_lib_LDPATH])dnl
+AC_SUBST([BOOST_LDPATH], [$Boost_lib_LDPATH])dnl
+AC_SUBST(AS_TR_CPP([BOOST_$1_LIBS]), [$Boost_lib_LIBS])dnl
+CPPFLAGS=$boost_save_CPPFLAGS
+AS_VAR_POPDEF([Boost_lib])dnl
+AS_VAR_POPDEF([Boost_lib_LDFLAGS])dnl
+AS_VAR_POPDEF([Boost_lib_LDPATH])dnl
+AS_VAR_POPDEF([Boost_lib_LIBS])dnl
+AC_LANG_POP([C++])dnl
+fi
+])
+
+
+# BOOST_FIND_LIB([LIB-NAME],
+# [PREFERRED-RT-OPT], [HEADER-NAME], [CXX-TEST],
+# [CXX-PROLOGUE])
+# --------------------------------------------------------------
+# Backward compatibility wrapper for BOOST_FIND_LIBS.
+AC_DEFUN([BOOST_FIND_LIB],
+[BOOST_FIND_LIBS([$1], $@)])
+
+
+# _BOOST_FIND_LIBS([LIB-NAME], [CANDIDATE-LIB-NAMES],
+# [PREFERRED-RT-OPT], [HEADER-NAME], [CXX-TEST],
+# [CXX-PROLOGUE])
+# --------------------------------------------------------------
+# Real implementation of BOOST_FIND_LIBS: rely on these local macros:
+# Boost_lib, Boost_lib_LDFLAGS, Boost_lib_LDPATH, Boost_lib_LIBS
+#
+# The algorithm is as follows: first look for a given library name
+# according to the user's PREFERRED-RT-OPT. For each library name, we
+# prefer to use the ones that carry the tag (toolset name). Each
+# library is searched through the various standard paths were Boost is
+# usually installed. If we can't find the standard variants, we try
+# to enforce -mt (for instance on MacOSX, libboost_thread.dylib
+# doesn't exist but there's -obviously- libboost_thread-mt.dylib).
+AC_DEFUN([_BOOST_FIND_LIBS],
+[Boost_lib=no
+ case "$3" in #(
+ (mt | mt-) boost_mt=-mt; boost_rtopt=;; #(
+ (mt* | mt-*) boost_mt=-mt; boost_rtopt=`expr "X$3" : 'Xmt-*\(.*\)'`;; #(
+ (*) boost_mt=; boost_rtopt=$3;;
esac
if test $enable_static_boost = yes; then
boost_rtopt="s$boost_rtopt"
fi
# Find the proper debug variant depending on what we've been asked to find.
case $boost_rtopt in #(
- *d*) boost_rt_d=$boost_rtopt;; #(
- *[[sgpn]]*) # Insert the `d' at the right place (in between `sg' and `pn')
+ (*d*) boost_rt_d=$boost_rtopt;; #(
+ (*[[sgpn]]*) # Insert the `d' at the right place (in between `sg' and `pn')
boost_rt_d=`echo "$boost_rtopt" | sed 's/\(s*g*\)\(p*n*\)/\1\2/'`;; #(
- *) boost_rt_d='-d';;
+ (*) boost_rt_d='-d';;
esac
# If the PREFERRED-RT-OPT are not empty, prepend a `-'.
test -n "$boost_rtopt" && boost_rtopt="-$boost_rtopt"
@@ -335,8 +395,8 @@ AC_CACHE_CHECK([for the Boost $1 library], [Boost_lib],
AC_MSG_ERROR([the libext variable is empty, did you invoke Libtool?])
boost_save_ac_objext=$ac_objext
# Generate the test file.
- AC_LANG_CONFTEST([AC_LANG_PROGRAM([#include <$3>
-$5], [$4])])
+ AC_LANG_CONFTEST([AC_LANG_PROGRAM([#include <$4>
+$6], [$5])])
dnl Optimization hacks: compiling C++ is slow, especially with Boost. What
dnl we're trying to do here is guess the right combination of link flags
dnl (LIBS / LDFLAGS) to use a given library. This can take several
@@ -358,21 +418,22 @@ dnl start the for loops).
[AC_MSG_ERROR([cannot compile a test that uses Boost $1])])
ac_objext=$boost_save_ac_objext
boost_failed_libs=
-# Don't bother to ident the 6 nested for loops, only the 2 innermost ones
-# matter.
+# Don't bother to ident the following nested for loops, only the 2
+# innermost ones matter.
+for boost_lib_ in $2; do
for boost_tag_ in -$boost_cv_lib_tag ''; do
for boost_ver_ in -$boost_cv_lib_version ''; do
for boost_mt_ in $boost_mt -mt ''; do
for boost_rtopt_ in $boost_rtopt '' -d; do
for boost_lib in \
- boost_$1$boost_tag_$boost_mt_$boost_rtopt_$boost_ver_ \
- boost_$1$boost_tag_$boost_rtopt_$boost_ver_ \
- boost_$1$boost_tag_$boost_mt_$boost_ver_ \
- boost_$1$boost_tag_$boost_ver_
+ boost_$boost_lib_$boost_tag_$boost_mt_$boost_rtopt_$boost_ver_ \
+ boost_$boost_lib_$boost_tag_$boost_rtopt_$boost_ver_ \
+ boost_$boost_lib_$boost_tag_$boost_mt_$boost_ver_ \
+ boost_$boost_lib_$boost_tag_$boost_ver_
do
# Avoid testing twice the same lib
case $boost_failed_libs in #(
- *@$boost_lib@*) continue;;
+ (*@$boost_lib@*) continue;;
esac
# If with_boost is empty, we'll search in /lib first, which is not quite
# right so instead we'll try to a location based on where the headers are.
@@ -382,14 +443,17 @@ for boost_rtopt_ in $boost_rtopt '' -d; do
/opt/local/lib* /usr/local/lib* /opt/lib* /usr/lib* \
"$with_boost" C:/Boost/lib /lib*
do
- test -e "$boost_ldpath" || continue
+ # Don't waste time with directories that don't exist.
+ if test x"$boost_ldpath" != x && test ! -e "$boost_ldpath"; then
+ continue
+ fi
boost_save_LDFLAGS=$LDFLAGS
# Are we looking for a static library?
case $boost_ldpath:$boost_rtopt_ in #(
- *?*:*s*) # Yes (Non empty boost_ldpath + s in rt opt)
+ (*?*:*s*) # Yes (Non empty boost_ldpath + s in rt opt)
Boost_lib_LIBS="$boost_ldpath/lib$boost_lib.$libext"
test -e "$Boost_lib_LIBS" || continue;; #(
- *) # No: use -lboost_foo to find the shared library.
+ (*) # No: use -lboost_foo to find the shared library.
Boost_lib_LIBS="-l$boost_lib";;
esac
boost_save_LIBS=$LIBS
@@ -403,9 +467,35 @@ dnl generated only once above (before we start the for loops).
LDFLAGS=$boost_save_LDFLAGS
LIBS=$boost_save_LIBS
if test x"$Boost_lib" = xyes; then
- Boost_lib_LDFLAGS="-L$boost_ldpath -Wl,-R$boost_ldpath"
+ # Check or used cached result of whether or not using -R or
+ # -rpath makes sense. Some implementations of ld, such as for
+ # Mac OSX, require -rpath but -R is the flag known to work on
+ # other systems. https://github.com/tsuna/boost.m4/issues/19
+ AC_CACHE_VAL([boost_cv_rpath_link_ldflag],
+ [case $boost_ldpath in
+ '') # Nothing to do.
+ boost_cv_rpath_link_ldflag=
+ boost_rpath_link_ldflag_found=yes;;
+ *)
+ for boost_cv_rpath_link_ldflag in -Wl,-R, -Wl,-rpath,; do
+ LDFLAGS="$boost_save_LDFLAGS -L$boost_ldpath $boost_cv_rpath_link_ldflag$boost_ldpath"
+ LIBS="$boost_save_LIBS $Boost_lib_LIBS"
+ _BOOST_AC_LINK_IFELSE([],
+ [boost_rpath_link_ldflag_found=yes
+ break],
+ [boost_rpath_link_ldflag_found=no])
+ done
+ ;;
+ esac
+ AS_IF([test "x$boost_rpath_link_ldflag_found" != "xyes"],
+ [AC_MSG_ERROR([Unable to determine whether to use -R or -rpath])])
+ LDFLAGS=$boost_save_LDFLAGS
+ LIBS=$boost_save_LIBS
+ ])
+ test x"$boost_ldpath" != x &&
+ Boost_lib_LDFLAGS="-L$boost_ldpath $boost_cv_rpath_link_ldflag$boost_ldpath"
Boost_lib_LDPATH="$boost_ldpath"
- break 6
+ break 7
else
boost_failed_libs="$boost_failed_libs@$boost_lib@"
fi
@@ -415,25 +505,10 @@ done
done
done
done
+done # boost_lib_
rm -f conftest.$ac_objext
])
-case $Boost_lib in #(
- no) _AC_MSG_LOG_CONFTEST
- AC_MSG_ERROR([cannot find the flags to link with Boost $1])
- ;;
-esac
-AC_SUBST(AS_TR_CPP([BOOST_$1_LDFLAGS]), [$Boost_lib_LDFLAGS])dnl
-AC_SUBST(AS_TR_CPP([BOOST_$1_LDPATH]), [$Boost_lib_LDPATH])dnl
-AC_SUBST([BOOST_LDPATH], [$Boost_lib_LDPATH])dnl
-AC_SUBST(AS_TR_CPP([BOOST_$1_LIBS]), [$Boost_lib_LIBS])dnl
-CPPFLAGS=$boost_save_CPPFLAGS
-AS_VAR_POPDEF([Boost_lib])dnl
-AS_VAR_POPDEF([Boost_lib_LDFLAGS])dnl
-AS_VAR_POPDEF([Boost_lib_LDPATH])dnl
-AS_VAR_POPDEF([Boost_lib_LIBS])dnl
-AC_LANG_POP([C++])dnl
-fi
-])# BOOST_FIND_LIB
+
# --------------------------------------- #
@@ -475,20 +550,20 @@ BOOST_FIND_HEADER([boost/asio.hpp])])
# BOOST_BIND()
# ------------
-# Look for Boost.Bind
+# Look for Boost.Bind.
BOOST_DEFUN([Bind],
[BOOST_FIND_HEADER([boost/bind.hpp])])
# BOOST_CHRONO()
-# ------------------
-# Look for Boost.Chrono
+# --------------
+# Look for Boost.Chrono.
BOOST_DEFUN([Chrono],
[# Do we have to check for Boost.System? This link-time dependency was
# added as of 1.35.0. If we have a version <1.35, we must not attempt to
# find Boost.System as it didn't exist by then.
if test $boost_major_version -ge 135; then
-BOOST_SYSTEM([$1])
+ BOOST_SYSTEM([$1])
fi # end of the Boost.System check.
boost_filesystem_save_LIBS=$LIBS
boost_filesystem_save_LDFLAGS=$LDFLAGS
@@ -499,14 +574,40 @@ BOOST_FIND_LIB([chrono], [$1],
[boost/chrono.hpp],
[boost::chrono::thread_clock d;])
if test $enable_static_boost = yes && test $boost_major_version -ge 135; then
- AC_SUBST([BOOST_FILESYSTEM_LIBS], ["$BOOST_FILESYSTEM_LIBS $BOOST_SYSTEM_LIBS"])
+ BOOST_CHRONO_LIBS="$BOOST_CHRONO_LIBS $BOOST_SYSTEM_LIBS"
fi
LIBS=$boost_filesystem_save_LIBS
LDFLAGS=$boost_filesystem_save_LDFLAGS
-
])# BOOST_CHRONO
+# BOOST_CONTEXT([PREFERRED-RT-OPT])
+# -----------------------------------
+# Look for Boost.Context. For the documentation of PREFERRED-RT-OPT, see the
+# documentation of BOOST_FIND_LIB above. This library was introduced in Boost
+# 1.51.0
+BOOST_DEFUN([Context],
+[BOOST_FIND_LIB([context], [$1],
+ [boost/context/all.hpp],[[
+// creates a stack
+void * stack_pointer = new void*[4096];
+std::size_t const size = sizeof(void*[4096]);
+
+// context fc uses f() as context function
+// fcontext_t is placed on top of context stack
+// a pointer to fcontext_t is returned
+fc = ctx::make_fcontext(stack_pointer, size, f);
+return ctx::jump_fcontext(&fcm, fc, 3) == 6;]],[dnl
+namespace ctx = boost::context;
+// context
+static ctx::fcontext_t fcm, *fc;
+// context-function
+static void f(intptr_t i) {
+ ctx::jump_fcontext(fc, &fcm, i * 2);
+}])
+])# BOOST_CONTEXT
+
+
# BOOST_CONVERSION()
# ------------------
# Look for Boost.Conversion (cast / lexical_cast)
@@ -516,6 +617,52 @@ BOOST_FIND_HEADER([boost/lexical_cast.hpp])
])# BOOST_CONVERSION
+# BOOST_COROUTINE([PREFERRED-RT-OPT])
+# -----------------------------------
+# Look for Boost.Coroutine. For the documentation of PREFERRED-RT-OPT, see the
+# documentation of BOOST_FIND_LIB above. This library was introduced in Boost
+# 1.53.0
+BOOST_DEFUN([Coroutine],
+[
+boost_coroutine_save_LIBS=$LIBS
+boost_coroutine_save_LDFLAGS=$LDFLAGS
+# Link-time dependency from coroutine to context
+BOOST_CONTEXT([$1])
+# Starting from Boost 1.55 a dependency on Boost.System is added
+if test $boost_major_version -ge 155; then
+ BOOST_SYSTEM([$1])
+fi
+m4_pattern_allow([^BOOST_(CONTEXT|SYSTEM)_(LIBS|LDFLAGS)])
+LIBS="$LIBS $BOOST_CONTEXT_LIBS $BOOST_SYSTEM_LIBS"
+LDFLAGS="$LDFLAGS $BOOST_CONTEXT_LDFLAGS"
+
+BOOST_FIND_LIB([coroutine], [$1],
+ [boost/coroutine/coroutine.hpp],
+ [boost::coroutines::coroutine< int(int) > coro; coro.empty();])
+
+# Link-time dependency from coroutine to context, existed only in 1.53, in 1.54
+# coroutine doesn't use context from its headers but from its library.
+if test $boost_major_version -eq 153 || test $enable_static_boost = yes && test $boost_major_version -ge 154; then
+ BOOST_COROUTINE_LIBS="$BOOST_COROUTINE_LIBS $BOOST_CONTEXT_LIBS"
+ BOOST_COROUTINE_LDFLAGS="$BOOST_COROUTINE_LDFLAGS $BOOST_CONTEXT_LDFLAGS"
+fi
+if test $enable_static_boost = yes && test $boost_major_version -ge 155; then
+ BOOST_COROUTINE_LIBS="$BOOST_COROUTINE_LIBS $BOOST_SYSTEM_LIBS"
+ BOOST_COROUTINE_LDFLAGS="$BOOST_COROUTINE_LDFLAGS $BOOST_SYSTEM_LDFLAGS"
+fi
+LIBS=$boost_coroutine_save_LIBS
+LDFLAGS=$boost_coroutine_save_LDFLAGS
+])# BOOST_COROUTINE
+
+
+# BOOST_CRC()
+# -----------
+# Look for Boost.CRC
+BOOST_DEFUN([CRC],
+[BOOST_FIND_HEADER([boost/crc.hpp])
+])# BOOST_CRC
+
+
# BOOST_DATE_TIME([PREFERRED-RT-OPT])
# -----------------------------------
# Look for Boost.Date_Time. For the documentation of PREFERRED-RT-OPT, see the
@@ -538,7 +685,7 @@ BOOST_DEFUN([Filesystem],
# added as of 1.35.0. If we have a version <1.35, we must not attempt to
# find Boost.System as it didn't exist by then.
if test $boost_major_version -ge 135; then
-BOOST_SYSTEM([$1])
+ BOOST_SYSTEM([$1])
fi # end of the Boost.System check.
boost_filesystem_save_LIBS=$LIBS
boost_filesystem_save_LDFLAGS=$LDFLAGS
@@ -548,23 +695,34 @@ LDFLAGS="$LDFLAGS $BOOST_SYSTEM_LDFLAGS"
BOOST_FIND_LIB([filesystem], [$1],
[boost/filesystem/path.hpp], [boost::filesystem::path p;])
if test $enable_static_boost = yes && test $boost_major_version -ge 135; then
- AC_SUBST([BOOST_FILESYSTEM_LIBS], ["$BOOST_FILESYSTEM_LIBS $BOOST_SYSTEM_LIBS"])
+ BOOST_FILESYSTEM_LIBS="$BOOST_FILESYSTEM_LIBS $BOOST_SYSTEM_LIBS"
fi
LIBS=$boost_filesystem_save_LIBS
LDFLAGS=$boost_filesystem_save_LDFLAGS
])# BOOST_FILESYSTEM
+# BOOST_FLYWEIGHT()
+# -----------------
+# Look for Boost.Flyweight.
+BOOST_DEFUN([Flyweight],
+[dnl There's a hidden dependency on pthreads.
+AC_REQUIRE([_BOOST_PTHREAD_FLAG])dnl
+BOOST_FIND_HEADER([boost/flyweight.hpp])
+AC_SUBST([BOOST_FLYWEIGHT_LIBS], [$boost_cv_pthread_flag])
+])
+
+
# BOOST_FOREACH()
# ---------------
-# Look for Boost.Foreach
+# Look for Boost.Foreach.
BOOST_DEFUN([Foreach],
[BOOST_FIND_HEADER([boost/foreach.hpp])])
# BOOST_FORMAT()
# --------------
-# Look for Boost.Format
+# Look for Boost.Format.
# Note: we can't check for boost/format/format_fwd.hpp because the header isn't
# standalone. It can't be compiled because it triggers the following error:
# boost/format/detail/config_macros.hpp:88: error: 'locale' in namespace 'std'
@@ -580,6 +738,14 @@ BOOST_DEFUN([Function],
[BOOST_FIND_HEADER([boost/function.hpp])])
+# BOOST_GEOMETRY()
+# ----------------
+# Look for Boost.Geometry (new since 1.47.0).
+BOOST_DEFUN([Geometry],
+[BOOST_FIND_HEADER([boost/geometry.hpp])
+])# BOOST_GEOMETRY
+
+
# BOOST_GRAPH([PREFERRED-RT-OPT])
# -------------------------------
# Look for Boost.Graphs. For the documentation of PREFERRED-RT-OPT, see the
@@ -615,9 +781,18 @@ BOOST_DEFUN([Lambda],
[BOOST_FIND_HEADER([boost/lambda/lambda.hpp])])
+# BOOST_LOCALE()
+# --------------
+# Look for Boost.Locale
+BOOST_DEFUN([Locale],
+[BOOST_FIND_LIB([locale], [$1],
+ [boost/locale.hpp],
+ [[boost::locale::generator gen; std::locale::global(gen(""));]])
+])# BOOST_LOCALE
+
# BOOST_LOG([PREFERRED-RT-OPT])
# -----------------------------
-# Look for Boost.Log For the documentation of PREFERRED-RT-OPT, see the
+# Look for Boost.Log. For the documentation of PREFERRED-RT-OPT, see the
# documentation of BOOST_FIND_LIB above.
BOOST_DEFUN([Log],
[BOOST_FIND_LIB([log], [$1],
@@ -628,12 +803,12 @@ BOOST_DEFUN([Log],
# BOOST_LOG_SETUP([PREFERRED-RT-OPT])
# -----------------------------------
-# Look for Boost.Log For the documentation of PREFERRED-RT-OPT, see the
+# Look for Boost.Log. For the documentation of PREFERRED-RT-OPT, see the
# documentation of BOOST_FIND_LIB above.
BOOST_DEFUN([Log_Setup],
[AC_REQUIRE([BOOST_LOG])dnl
BOOST_FIND_LIB([log_setup], [$1],
- [boost/log/utility/init/from_settings.hpp],
+ [boost/log/utility/setup/from_settings.hpp],
[boost::log::basic_settings<char> bs; bs.empty();])
])# BOOST_LOG_SETUP
@@ -650,6 +825,29 @@ BOOST_DEFUN([Math],
[BOOST_FIND_HEADER([boost/math/special_functions.hpp])])
+# BOOST_MPI([PREFERRED-RT-OPT])
+# -------------------------------
+# Look for Boost MPI. For the documentation of PREFERRED-RT-OPT, see the
+# documentation of BOOST_FIND_LIB above. Uses MPICXX variable if it is
+# set, otherwise tries CXX
+#
+BOOST_DEFUN([MPI],
+[boost_save_CXX=${CXX}
+boost_save_CXXCPP=${CXXCPP}
+if test x"${MPICXX}" != x; then
+ CXX=${MPICXX}
+ CXXCPP="${MPICXX} -E"
+fi
+BOOST_FIND_LIB([mpi], [$1],
+ [boost/mpi.hpp],
+ [int argc = 0;
+ char **argv = 0;
+ boost::mpi::environment env(argc,argv);])
+CXX=${boost_save_CXX}
+CXXCPP=${boost_save_CXXCPP}
+])# BOOST_MPI
+
+
# BOOST_MULTIARRAY()
# ------------------
# Look for Boost.MultiArray
@@ -657,6 +855,14 @@ BOOST_DEFUN([MultiArray],
[BOOST_FIND_HEADER([boost/multi_array.hpp])])
+# BOOST_NUMERIC_UBLAS()
+# --------------------------
+# Look for Boost.NumericUblas (Basic Linear Algebra)
+BOOST_DEFUN([Numeric_Ublas],
+[BOOST_FIND_HEADER([boost/numeric/ublas/vector.hpp])
+])# BOOST_NUMERIC_UBLAS
+
+
# BOOST_NUMERIC_CONVERSION()
# --------------------------
# Look for Boost.NumericConversion (policy-based numeric conversion)
@@ -679,6 +885,12 @@ BOOST_DEFUN([Preprocessor],
[BOOST_FIND_HEADER([boost/preprocessor/repeat.hpp])])
+# BOOST_RANGE()
+# --------------------
+# Look for Boost.Range
+BOOST_DEFUN([Range],
+[BOOST_FIND_HEADER([boost/range/adaptors.hpp])])
+
# BOOST_UNORDERED()
# -----------------
# Look for Boost.Unordered
@@ -725,9 +937,9 @@ BOOST_DEFUN([Python],
_BOOST_PYTHON_CONFIG([LDFLAGS], [ldflags])
_BOOST_PYTHON_CONFIG([LIBS], [libs])
m4_pattern_allow([^BOOST_PYTHON_MODULE$])dnl
-BOOST_FIND_LIB([python], [$1],
- [boost/python.hpp],
- [], [BOOST_PYTHON_MODULE(empty) {}])
+BOOST_FIND_LIBS([python], [python python3], [$1],
+ [boost/python.hpp],
+ [], [BOOST_PYTHON_MODULE(empty) {}])
CPPFLAGS=$boost_python_save_CPPFLAGS
LDFLAGS=$boost_python_save_LDFLAGS
LIBS=$boost_python_save_LIBS
@@ -775,6 +987,14 @@ BOOST_DEFUN([Signals],
])# BOOST_SIGNALS
+# BOOST_SIGNALS2()
+# ----------------
+# Look for Boost.Signals2 (new since 1.39.0).
+BOOST_DEFUN([Signals2],
+[BOOST_FIND_HEADER([boost/signals2.hpp])
+])# BOOST_SIGNALS2
+
+
# BOOST_SMART_PTR()
# -----------------
# Look for Boost.SmartPtr
@@ -825,19 +1045,18 @@ BOOST_FIND_LIB([unit_test_framework], [$1],
])# BOOST_TEST
-# BOOST_THREADS([PREFERRED-RT-OPT])
+# BOOST_THREAD([PREFERRED-RT-OPT])
# ---------------------------------
# Look for Boost.Thread. For the documentation of PREFERRED-RT-OPT, see the
# documentation of BOOST_FIND_LIB above.
-# FIXME: Provide an alias "BOOST_THREAD".
-BOOST_DEFUN([Threads],
+BOOST_DEFUN([Thread],
[dnl Having the pthread flag is required at least on GCC3 where
dnl boost/thread.hpp would complain if we try to compile without
dnl -pthread on GNU/Linux.
AC_REQUIRE([_BOOST_PTHREAD_FLAG])dnl
-boost_threads_save_LIBS=$LIBS
-boost_threads_save_LDFLAGS=$LDFLAGS
-boost_threads_save_CPPFLAGS=$CPPFLAGS
+boost_thread_save_LIBS=$LIBS
+boost_thread_save_LDFLAGS=$LDFLAGS
+boost_thread_save_CPPFLAGS=$CPPFLAGS
# Link-time dependency from thread to system was added as of 1.49.0.
if test $boost_major_version -ge 149; then
BOOST_SYSTEM([$1])
@@ -845,36 +1064,26 @@ fi # end of the Boost.System check.
m4_pattern_allow([^BOOST_SYSTEM_(LIBS|LDFLAGS)$])dnl
LIBS="$LIBS $BOOST_SYSTEM_LIBS $boost_cv_pthread_flag"
LDFLAGS="$LDFLAGS $BOOST_SYSTEM_LDFLAGS"
-# Yes, we *need* to put the -pthread thing in CPPFLAGS because with GCC3,
-# boost/thread.hpp will trigger a #error if -pthread isn't used:
-# boost/config/requires_threads.hpp:47:5: #error "Compiler threading support
-# is not turned on. Please set the correct command line options for
-# threading: -pthread (Linux), -pthreads (Solaris) or -mthreads (Mingw32)"
CPPFLAGS="$CPPFLAGS $boost_cv_pthread_flag"
# When compiling for the Windows platform, the threads library is named
# differently.
case $host_os in
- (*mingw*)
- BOOST_FIND_LIB([thread_win32], [$1],
- [boost/thread.hpp], [boost::thread t; boost::mutex m;])
- BOOST_THREAD_LDFLAGS=$BOOST_THREAD_WIN32_LDFLAGS
- BOOST_THREAD_LDPATH=$BOOST_THREAD_WIN32_LDPATH
- BOOST_THREAD_LIBS=$BOOST_THREAD_WIN32_LIBS
- ;;
- (*)
- BOOST_FIND_LIB([thread], [$1],
- [boost/thread.hpp], [boost::thread t; boost::mutex m;])
- ;;
+ (*mingw*) boost_thread_lib_ext=_win32;;
esac
+BOOST_FIND_LIBS([thread], [thread$boost_thread_lib_ext],
+ [$1],
+ [boost/thread.hpp], [boost::thread t; boost::mutex m;])
BOOST_THREAD_LIBS="$BOOST_THREAD_LIBS $BOOST_SYSTEM_LIBS $boost_cv_pthread_flag"
BOOST_THREAD_LDFLAGS="$BOOST_SYSTEM_LDFLAGS"
BOOST_CPPFLAGS="$BOOST_CPPFLAGS $boost_cv_pthread_flag"
-LIBS=$boost_threads_save_LIBS
-LDFLAGS=$boost_threads_save_LDFLAGS
-CPPFLAGS=$boost_threads_save_CPPFLAGS
-])# BOOST_THREADS
+LIBS=$boost_thread_save_LIBS
+LDFLAGS=$boost_thread_save_LDFLAGS
+CPPFLAGS=$boost_thread_save_CPPFLAGS
+])# BOOST_THREAD
+
+AU_ALIAS([BOOST_THREADS], [BOOST_THREAD])
# BOOST_TOKENIZER()
@@ -923,10 +1132,23 @@ BOOST_DEFUN([Variant],
BOOST_FIND_HEADER([boost/variant.hpp])])
+# BOOST_POINTER_CONTAINER()
+# ------------------------
+# Look for Boost.PointerContainer
+BOOST_DEFUN([Pointer_Container],
+[BOOST_FIND_HEADER([boost/ptr_container/ptr_deque.hpp])
+BOOST_FIND_HEADER([boost/ptr_container/ptr_list.hpp])
+BOOST_FIND_HEADER([boost/ptr_container/ptr_vector.hpp])
+BOOST_FIND_HEADER([boost/ptr_container/ptr_array.hpp])
+BOOST_FIND_HEADER([boost/ptr_container/ptr_set.hpp])
+BOOST_FIND_HEADER([boost/ptr_container/ptr_map.hpp])
+])# BOOST_POINTER_CONTAINER
+
+
# BOOST_WAVE([PREFERRED-RT-OPT])
# ------------------------------
# NOTE: If you intend to use Wave/Spirit with thread support, make sure you
-# call BOOST_THREADS first.
+# call BOOST_THREAD first.
# Look for Boost.Wave. For the documentation of PREFERRED-RT-OPT, see the
# documentation of BOOST_FIND_LIB above.
BOOST_DEFUN([Wave],
@@ -961,8 +1183,16 @@ BOOST_DEFUN([Xpressive],
# _BOOST_PTHREAD_FLAG()
# ---------------------
-# Internal helper for BOOST_THREADS. Based on ACX_PTHREAD:
-# http://autoconf-archive.cryp.to/acx_pthread.html
+# Internal helper for BOOST_THREAD. Computes boost_cv_pthread_flag
+# which must be used in CPPFLAGS and LIBS.
+#
+# Yes, we *need* to put the -pthread thing in CPPFLAGS because with GCC3,
+# boost/thread.hpp will trigger a #error if -pthread isn't used:
+# boost/config/requires_threads.hpp:47:5: #error "Compiler threading support
+# is not turned on. Please set the correct command line options for
+# threading: -pthread (Linux), -pthreads (Solaris) or -mthreads (Mingw32)"
+#
+# Based on ACX_PTHREAD: http://autoconf-archive.cryp.to/acx_pthread.html
AC_DEFUN([_BOOST_PTHREAD_FLAG],
[AC_REQUIRE([AC_PROG_CXX])dnl
AC_REQUIRE([AC_CANONICAL_HOST])dnl
@@ -1030,6 +1260,14 @@ AC_LANG_POP([C++])dnl
m4_define([_BOOST_gcc_test],
["defined __GNUC__ && __GNUC__ == $1 && __GNUC_MINOR__ == $2 && !defined __ICC @ gcc$1$2"])dnl
+# _BOOST_mingw_test(MAJOR, MINOR)
+# -----------------------------
+# Internal helper for _BOOST_FIND_COMPILER_TAG.
+m4_define([_BOOST_mingw_test],
+["defined __GNUC__ && __GNUC__ == $1 && __GNUC_MINOR__ == $2 && !defined __ICC && \
+ (defined WIN32 || defined WINNT || defined _WIN32 || defined __WIN32 \
+ || defined __WIN32__ || defined __WINNT || defined __WINNT__) @ mgw$1$2"])dnl
+
# _BOOST_FIND_COMPILER_TAG()
# --------------------------
@@ -1039,7 +1277,8 @@ m4_define([_BOOST_gcc_test],
AC_DEFUN([_BOOST_FIND_COMPILER_TAG],
[AC_REQUIRE([AC_PROG_CXX])dnl
AC_REQUIRE([AC_CANONICAL_HOST])dnl
-AC_CACHE_CHECK([for the toolset name used by Boost for $CXX], [boost_cv_lib_tag],
+AC_CACHE_CHECK([for the toolset name used by Boost for $CXX],
+ [boost_cv_lib_tag],
[boost_cv_lib_tag=unknown
if test x$boost_cv_inc_path != xno; then
AC_LANG_PUSH([C++])dnl
@@ -1057,14 +1296,29 @@ if test x$boost_cv_inc_path != xno; then
# I'm not sure about my test for `il' (be careful: Intel's ICC pre-defines
# the same defines as GCC's).
for i in \
+ _BOOST_mingw_test(5, 0) \
+ _BOOST_gcc_test(5, 0) \
+ _BOOST_mingw_test(4, 10) \
+ _BOOST_gcc_test(4, 10) \
+ _BOOST_mingw_test(4, 9) \
+ _BOOST_gcc_test(4, 9) \
+ _BOOST_mingw_test(4, 8) \
_BOOST_gcc_test(4, 8) \
+ _BOOST_mingw_test(4, 7) \
_BOOST_gcc_test(4, 7) \
+ _BOOST_mingw_test(4, 6) \
_BOOST_gcc_test(4, 6) \
+ _BOOST_mingw_test(4, 5) \
_BOOST_gcc_test(4, 5) \
+ _BOOST_mingw_test(4, 4) \
_BOOST_gcc_test(4, 4) \
+ _BOOST_mingw_test(4, 3) \
_BOOST_gcc_test(4, 3) \
+ _BOOST_mingw_test(4, 2) \
_BOOST_gcc_test(4, 2) \
+ _BOOST_mingw_test(4, 1) \
_BOOST_gcc_test(4, 1) \
+ _BOOST_mingw_test(4, 0) \
_BOOST_gcc_test(4, 0) \
"defined __GNUC__ && __GNUC__ == 3 && !defined __ICC \
&& (defined WIN32 || defined WINNT || defined _WIN32 || defined __WIN32 \
@@ -1130,6 +1384,7 @@ fi])dnl end of AC_CACHE_CHECK
# Thread) flavors of Boost. Sets boost_guess_use_mt accordingly.
AC_DEFUN([_BOOST_GUESS_WHETHER_TO_USE_MT],
[# Check whether we do better use `mt' even though we weren't ask to.
+AC_LANG_PUSH([C++])dnl
AC_COMPILE_IFELSE([AC_LANG_PROGRAM([[
#if defined _REENTRANT || defined _MT || defined __MT__
/* use -mt */
@@ -1137,6 +1392,7 @@ AC_COMPILE_IFELSE([AC_LANG_PROGRAM([[
# error MT not needed
#endif
]])], [boost_guess_use_mt=:], [boost_guess_use_mt=false])
+AC_LANG_POP([C++])dnl
])
# _BOOST_AC_LINK_IFELSE(PROGRAM, [ACTION-IF-TRUE], [ACTION-IF-FALSE])
@@ -1160,11 +1416,11 @@ boost_use_source=:
test -f conftest.$ac_objext && ac_ext=$ac_objext && boost_use_source=false &&
_AS_ECHO_LOG([re-using the existing conftest.$ac_objext])
AS_IF([_AC_DO_STDERR($ac_link) && {
- test -z "$ac_[]_AC_LANG_ABBREV[]_werror_flag" ||
- test ! -s conftest.err
+ test -z "$ac_[]_AC_LANG_ABBREV[]_werror_flag" ||
+ test ! -s conftest.err
} && test -s conftest$ac_exeext && {
- test "$cross_compiling" = yes ||
- $as_executable_p conftest$ac_exeext
+ test "$cross_compiling" = yes ||
+ $as_executable_p conftest$ac_exeext
dnl FIXME: use AS_TEST_X instead when 2.61 is widespread enough.
}],
[$2],
diff --git a/src/Makefile.am b/src/Makefile.am
index cf3c6f8..b45aff3 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -28,7 +28,7 @@
-AM_CPPFLAGS = @GDALCFLAGS@ @BOOST_CPPFLAGS@ -Wall -Wextra -O2 -g
+AM_CPPFLAGS = @GDALCFLAGS@ @BOOST_CPPFLAGS@ -Wall -Wextra -O3 -g
LIBS = @GDALLIBS@
bin_PROGRAMS = gdal_raw2geotiff gdal_dem2rgb gdal_list_corners gdal_trace_outline gdal_contrast_stretch gdal_landsat_pansharp gdal_wkt_to_mask gdal_merge_simple gdal_merge_vrt gdal_get_projected_bounds gdal_make_ndv_mask
@@ -39,13 +39,13 @@ bin_PROGRAMS = gdal_raw2geotiff gdal_dem2rgb gdal_list_corners gdal_trace_outlin
gdal_raw2geotiff_SOURCES = gdal_raw2geotiff.cc common.cc
palette.o: default_palette.h
-gdal_dem2rgb_SOURCES = gdal_dem2rgb.cc common.cc georef.cc ndv.cc palette.cc
+gdal_dem2rgb_SOURCES = gdal_dem2rgb.cc common.cc georef.cc ndv.cc palette.cc datatype_conversion.cc
-gdal_list_corners_SOURCES = gdal_list_corners.cc common.cc polygon.cc polygon-rasterizer.cc debugplot.cc georef.cc mask.cc rectangle_finder.cc ndv.cc
+gdal_list_corners_SOURCES = gdal_list_corners.cc common.cc polygon.cc polygon-rasterizer.cc debugplot.cc georef.cc mask.cc rectangle_finder.cc ndv.cc datatype_conversion.cc
-gdal_trace_outline_SOURCES = gdal_trace_outline.cc common.cc polygon.cc polygon-rasterizer.cc debugplot.cc georef.cc mask.cc mask-tracer.cc beveler.cc dp.cc ndv.cc excursion_pincher2.cc
+gdal_trace_outline_SOURCES = gdal_trace_outline.cc common.cc polygon.cc polygon-rasterizer.cc debugplot.cc georef.cc mask.cc mask-tracer.cc beveler.cc dp.cc ndv.cc excursion_pincher2.cc raster_features.cc datatype_conversion.cc
-gdal_contrast_stretch_SOURCES = gdal_contrast_stretch.cc common.cc ndv.cc
+gdal_contrast_stretch_SOURCES = gdal_contrast_stretch.cc common.cc ndv.cc datatype_conversion.cc
gdal_landsat_pansharp_SOURCES = gdal_landsat_pansharp.cc common.cc
@@ -57,7 +57,7 @@ gdal_merge_simple_SOURCES = gdal_merge_simple.cc common.cc
gdal_merge_vrt_SOURCES = gdal_merge_vrt.cc common.cc
-gdal_make_ndv_mask_SOURCES = gdal_make_ndv_mask.cc common.cc ndv.cc mask.cc debugplot.cc
+gdal_make_ndv_mask_SOURCES = gdal_make_ndv_mask.cc common.cc ndv.cc mask.cc debugplot.cc datatype_conversion.cc
lint:
cpplint.py --filter=-whitespace,-readability/streams,-build/header_guard,-build/include_order,-readability/multiline_string \
diff --git a/src/common.cc b/src/common.cc
index 9b893ff..59b1c98 100644
--- a/src/common.cc
+++ b/src/common.cc
@@ -31,6 +31,8 @@ This code was developed by Dan Stahlke for the Geographic Information Network of
#include "common.h"
+namespace dangdal {
+
int VERBOSE = 0;
void fatal_error(const std::string &s) {
@@ -59,3 +61,5 @@ std::vector<std::string> argv_to_list(int argc, char **argv) {
}
return ret;
}
+
+} // namespace dangdal
diff --git a/src/common.h b/src/common.h
index f955eb2..59f7fb7 100644
--- a/src/common.h
+++ b/src/common.h
@@ -76,10 +76,14 @@ This code was developed by Dan Stahlke for the Geographic Information Network of
#define D2R (M_PI / 180.0)
+namespace dangdal {
+
extern int VERBOSE;
void fatal_error(const std::string &s) __attribute__((noreturn));
void fatal_error(const char *s, ...) __attribute__((noreturn, format(printf, 1, 2)));
std::vector<std::string> argv_to_list(int argc, char **argv);
+} // namespace dangdal
+
#endif // ifndef DANGDAL_COMMON_H
diff --git a/src/datatype_conversion.cc b/src/datatype_conversion.cc
new file mode 100644
index 0000000..11a2b09
--- /dev/null
+++ b/src/datatype_conversion.cc
@@ -0,0 +1,65 @@
+#include <gdal.h>
+
+#include "datatype_conversion.h"
+
+namespace dangdal {
+
+int32_t gdal_scalar_to_int32(void *p, GDALDataType dt) {
+ int32_t ret;
+ GDALCopyWords(p, dt, 0, &ret, GDT_Int32, 0, 1);
+ return ret;
+}
+
+double gdal_scalar_to_double(void *p, GDALDataType dt) {
+ double ret;
+ GDALCopyWords(p, dt, 0, &ret, GDT_Float64, 0, 1);
+ return ret;
+}
+
+template <typename T>
+struct NaN_Checker {
+ static bool isnan(const void *) { return false; }
+};
+
+template <>
+struct NaN_Checker<double> {
+ static bool isnan(const void *p) {
+ return std::isnan(*reinterpret_cast<const double *>(p));
+ }
+};
+
+template <>
+struct NaN_Checker<float> {
+ static bool isnan(const void *p) {
+ return std::isnan(*reinterpret_cast<const float *>(p));
+ }
+};
+
+template <>
+struct NaN_Checker<std::complex<double> > {
+ static bool isnan(const void *p) {
+ return
+ std::isnan(reinterpret_cast<const std::complex<double> *>(p)->real()) ||
+ std::isnan(reinterpret_cast<const std::complex<double> *>(p)->imag());
+ }
+};
+
+template <>
+struct NaN_Checker<std::complex<float> > {
+ static bool isnan(const void *p) {
+ return
+ std::isnan(reinterpret_cast<const std::complex<float> *>(p)->real()) ||
+ std::isnan(reinterpret_cast<const std::complex<float> *>(p)->imag());
+ }
+};
+
+template <typename T>
+static bool nan_check(const void *p) {
+ return NaN_Checker<T>::isnan(p);
+}
+
+bool gdal_scalar_pointer_isnan(const void *p, GDALDataType dt) {
+ return DANGDAL_RUNTIME_TEMPLATE(dt, nan_check, p);
+}
+
+} // namespace dangdal
diff --git a/src/datatype_conversion.h b/src/datatype_conversion.h
new file mode 100644
index 0000000..77126ac
--- /dev/null
+++ b/src/datatype_conversion.h
@@ -0,0 +1,49 @@
+#ifndef DANGDAL_DATATYPE_CONVERSION_H
+#define DANGDAL_DATATYPE_CONVERSION_H
+
+#include <stdint.h>
+
+#include <complex>
+#include <stdexcept>
+
+namespace dangdal {
+
+// This lets you convert a GDALDataType to a template argument.
+#define DANGDAL_RUNTIME_TEMPLATE(dt, fn, ...) ( \
+ (dt == GDT_Byte ) ? fn< uint8_t>(__VA_ARGS__) : \
+ (dt == GDT_UInt16 ) ? fn<uint16_t>(__VA_ARGS__) : \
+ (dt == GDT_Int16 ) ? fn< int16_t>(__VA_ARGS__) : \
+ (dt == GDT_UInt32 ) ? fn<uint32_t>(__VA_ARGS__) : \
+ (dt == GDT_Int32 ) ? fn< int32_t>(__VA_ARGS__) : \
+ (dt == GDT_Float32 ) ? fn< float>(__VA_ARGS__) : \
+ (dt == GDT_Float64 ) ? fn< double>(__VA_ARGS__) : \
+ (dt == GDT_CInt16 ) ? fn<std::complex<int16_t> >(__VA_ARGS__) : \
+ (dt == GDT_CInt32 ) ? fn<std::complex<int32_t> >(__VA_ARGS__) : \
+ (dt == GDT_CFloat32) ? fn<std::complex< float> >(__VA_ARGS__) : \
+ (dt == GDT_CFloat64) ? fn<std::complex< double> >(__VA_ARGS__) : \
+ throw(std::invalid_argument("unrecognized datatype")) \
+)
+
+template <typename T>
+struct GetGDALDataTypeFor { static const GDALDataType t = GDT_Unknown; };
+template <> struct GetGDALDataTypeFor< uint8_t> { static const GDALDataType t = GDT_Byte ; };
+template <> struct GetGDALDataTypeFor<uint16_t> { static const GDALDataType t = GDT_UInt16 ; };
+template <> struct GetGDALDataTypeFor< int16_t> { static const GDALDataType t = GDT_Int16 ; };
+template <> struct GetGDALDataTypeFor<uint32_t> { static const GDALDataType t = GDT_UInt32 ; };
+template <> struct GetGDALDataTypeFor< int32_t> { static const GDALDataType t = GDT_Int32 ; };
+template <> struct GetGDALDataTypeFor< float> { static const GDALDataType t = GDT_Float32; };
+template <> struct GetGDALDataTypeFor< double> { static const GDALDataType t = GDT_Float64; };
+template <> struct GetGDALDataTypeFor<std::complex<int16_t> > { static const GDALDataType t = GDT_CInt16 ; };
+template <> struct GetGDALDataTypeFor<std::complex<int32_t> > { static const GDALDataType t = GDT_CInt32 ; };
+template <> struct GetGDALDataTypeFor<std::complex< float> > { static const GDALDataType t = GDT_CFloat32; };
+template <> struct GetGDALDataTypeFor<std::complex< double> > { static const GDALDataType t = GDT_CFloat64; };
+
+// unfortunately, GDALCopyWords doesn't allow const input
+int32_t gdal_scalar_to_int32(void *p, GDALDataType dt);
+double gdal_scalar_to_double(void *p, GDALDataType dt);
+
+bool gdal_scalar_pointer_isnan(const void *p, GDALDataType dt);
+
+} // namespace dangdal
+
+#endif // DANGDAL_DATATYPE_CONVERSION_H
diff --git a/src/default_palette.h b/src/default_palette.h
index a7db1e1..e7caff9 100644
--- a/src/default_palette.h
+++ b/src/default_palette.h
@@ -1,3 +1,34 @@
+/*
+Copyright (c) 2013, Regents of the University of Alaska
+
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
+
+ * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
+ * Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
+ * Neither the name of the Geographic Information Network of Alaska nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
+CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+This code was developed by Kevin Engle for the Geographic Information Network of Alaska.
+*/
+
+
+
+namespace dangdal {
+
+// Palette by Kevin Engle.
const char *DEFAULT_PALETTE[] = {
"NaN 0 0 0",
"-10200.000000 0 0 255",
@@ -163,3 +194,5 @@ const char *DEFAULT_PALETTE[] = {
"6900.000000 104 46 3",
"7000.000000 104 46 3",
NULL };
+
+} // namespace dangdal
diff --git a/src/dp.cc b/src/dp.cc
index 6a85627..ddb19f5 100644
--- a/src/dp.cc
+++ b/src/dp.cc
@@ -27,6 +27,8 @@ This code was developed by Dan Stahlke for the Geographic Information Network of
#include <vector>
+#include <utility>
+#include <cassert>
#include "common.h"
#include "polygon.h"
@@ -333,11 +335,33 @@ static inline int segs_cross(
1);
}
+// This allows rapidly finding which segments intersect each other.
+// The std::pair consists of the ring and segment indices.
+BboxBinarySpacePartition<std::pair<size_t, size_t> > get_bsp_for_reduced_rings(
+ const Mpoly &mpoly, std::vector<ReducedRing> &reduced_rings
+) {
+ typedef std::pair<size_t, size_t> segptr_t;
+ std::vector<std::pair<Bbox, segptr_t> > items;
+
+ for(size_t ring_idx=0; ring_idx < mpoly.rings.size(); ring_idx++) {
+ const Ring &ring = mpoly.rings[ring_idx];
+ const ReducedRing &reduced = reduced_rings[ring_idx];
+ for(size_t seg_idx=0; seg_idx < reduced.segs.size(); seg_idx++) {
+ Bbox bbox = reduced.segs[seg_idx].get_bbox(ring);
+ items.push_back(std::make_pair(bbox, segptr_t(ring_idx, seg_idx)));
+ }
+ }
+
+ return BboxBinarySpacePartition<segptr_t>(items);
+}
+
void fix_topology(const Mpoly &mpoly, std::vector<ReducedRing> &reduced_rings) {
const double firsthalf_progress = 0.5;
printf("Fixing topology: ");
fflush(stdout);
+ assert(mpoly.rings.size() == reduced_rings.size());
+
// initialize problem arrays
std::vector<std::vector<bool> > mp_problems(mpoly.rings.size());
for(size_t r1_idx=0; r1_idx < mpoly.rings.size(); r1_idx++) {
@@ -345,31 +369,34 @@ void fix_topology(const Mpoly &mpoly, std::vector<ReducedRing> &reduced_rings) {
mp_problems[r1_idx].resize(rring.segs.size(), 0);
}
- std::vector<Bbox> bboxes = mpoly.getRingBboxes();
+ int num_problems = 0;
+ {
+ const BboxBinarySpacePartition<std::pair<size_t, size_t> > bsp =
+ get_bsp_for_reduced_rings(mpoly, reduced_rings);
- // flag segments that cross
- int have_problems = 0;
- for(size_t r1_idx=0; r1_idx < mpoly.rings.size(); r1_idx++) {
- GDALTermProgress(firsthalf_progress*
- pow((double)r1_idx / (double)mpoly.rings.size(), 2), NULL, NULL);
- const Ring &c1 = mpoly.rings[r1_idx];
- const ReducedRing &r1 = reduced_rings[r1_idx];
- std::vector<bool> &p1 = mp_problems[r1_idx];
- const Bbox bbox1 = bboxes[r1_idx];
- if(bbox1.empty) continue;
- for(size_t r2_idx=0; r2_idx < mpoly.rings.size(); r2_idx++) {
- if(r2_idx > r1_idx) continue; // symmetry optimization
-
- const Bbox bbox2 = bboxes[r2_idx];
- if(is_disjoint(bbox1, bbox2)) continue;
-
- const Ring &c2 = mpoly.rings[r2_idx];
- const ReducedRing &r2 = reduced_rings[r2_idx];
- std::vector<bool> &p2 = mp_problems[r2_idx];
+ // flag segments that cross
+ for(size_t r1_idx=0; r1_idx < mpoly.rings.size(); r1_idx++) {
+ GDALTermProgress(firsthalf_progress*
+ pow((double)r1_idx / (double)mpoly.rings.size(), 2), NULL, NULL);
+ const Ring &c1 = mpoly.rings[r1_idx];
+ const ReducedRing &r1 = reduced_rings[r1_idx];
+ std::vector<bool> &p1 = mp_problems[r1_idx];
for(size_t seg1_idx=0; seg1_idx < r1.segs.size(); seg1_idx++) {
- for(size_t seg2_idx=0; seg2_idx < r2.segs.size(); seg2_idx++) {
- if(r2_idx == r1_idx && seg2_idx > seg1_idx) continue; // symmetry optimization
+ Bbox seg1_bbox = r1.segs[seg1_idx].get_bbox(c1);
+ std::vector<std::pair<size_t, size_t> > intersecting_segments =
+ bsp.get_intersecting_items(seg1_bbox);
+ for(size_t i_s_idx=0; i_s_idx < intersecting_segments.size(); i_s_idx++) {
+ size_t r2_idx = intersecting_segments[i_s_idx].first;
+ size_t seg2_idx = intersecting_segments[i_s_idx].second;
+
+ // symmetry optimization
+ if(r2_idx > r1_idx) continue;
+ if(r2_idx == r1_idx && seg2_idx > seg1_idx) continue;
+
+ const Ring &c2 = mpoly.rings[r2_idx];
+ const ReducedRing &r2 = reduced_rings[r2_idx];
+ std::vector<bool> &p2 = mp_problems[r2_idx];
int crosses = segs_cross(r1_idx==r2_idx,
c1, r1.segs[seg1_idx], c2, r2.segs[seg2_idx]);
@@ -378,23 +405,23 @@ void fix_topology(const Mpoly &mpoly, std::vector<ReducedRing> &reduced_rings) {
// r1_idx, seg1_idx, r2_idx, seg2_idx);
p1[seg1_idx] = 1;
p2[seg2_idx] = 1;
- have_problems += 2;
+ num_problems += 2;
}
- } // seg loop
- } // seg loop
- } // ring loop
- } // ring loop
+ } // ring2/seg2 loop
+ } // seg1 loop
+ } // ring1 loop
+ }
double progress = firsthalf_progress;
GDALTermProgress(progress, NULL, NULL);
- if(have_problems) {
- if(VERBOSE) printf("fixing %d crossed segments from reduction\n", have_problems/2);
+ if(num_problems) {
+ if(VERBOSE) printf("fixing %d crossed segments from reduction\n", num_problems/2);
}
int did_something = 1;
- while(have_problems && did_something) {
- //printf("%d crossings to fix\n", have_problems/2);
+ while(num_problems && did_something) {
+ //printf("%d crossings to fix\n", num_problems/2);
did_something = 0;
// subdivide problem segments
for(size_t r1_idx=0; r1_idx < mpoly.rings.size(); r1_idx++) {
@@ -417,7 +444,10 @@ void fix_topology(const Mpoly &mpoly, std::vector<ReducedRing> &reduced_rings) {
} // seg loop
} // ring loop
- have_problems = 0;
+ const BboxBinarySpacePartition<std::pair<size_t, size_t> > bsp =
+ get_bsp_for_reduced_rings(mpoly, reduced_rings);
+
+ num_problems = 0;
// now test for resolved problems and new problems
for(size_t r1_idx=0; r1_idx < mpoly.rings.size(); r1_idx++) {
{
@@ -428,48 +458,51 @@ void fix_topology(const Mpoly &mpoly, std::vector<ReducedRing> &reduced_rings) {
const Ring &c1 = mpoly.rings[r1_idx];
const ReducedRing &r1 = reduced_rings[r1_idx];
- const Bbox bbox1 = bboxes[r1_idx];
std::vector<bool> &p1 = mp_problems[r1_idx];
for(size_t seg1_idx=0; seg1_idx < r1.segs.size(); seg1_idx++) {
if(!p1[seg1_idx]) continue;
p1[seg1_idx] = 0;
- for(size_t r2_idx=0; r2_idx < mpoly.rings.size(); r2_idx++) {
+
+ Bbox seg1_bbox = r1.segs[seg1_idx].get_bbox(c1);
+ std::vector<std::pair<size_t, size_t> > intersecting_segments =
+ bsp.get_intersecting_items(seg1_bbox);
+ for(size_t i_s_idx=0; i_s_idx < intersecting_segments.size(); i_s_idx++) {
+ size_t r2_idx = intersecting_segments[i_s_idx].first;
+ size_t seg2_idx = intersecting_segments[i_s_idx].second;
+
const Ring &c2 = mpoly.rings[r2_idx];
const ReducedRing &r2 = reduced_rings[r2_idx];
- const Bbox bbox2 = bboxes[r2_idx];
- if(is_disjoint(bbox1, bbox2)) continue;
- for(size_t seg2_idx=0; seg2_idx < r2.segs.size(); seg2_idx++) {
- int crosses = segs_cross(r1_idx==r2_idx,
- c1, r1.segs[seg1_idx], c2, r2.segs[seg2_idx]);
- if(crosses) {
- if(VERBOSE) {
- printf("found a crossing (still): %zd,%zd,%zd,%zd (%f,%f)-(%f,%f) (%f,%f)-(%f,%f)\n",
- r1_idx, seg1_idx, r2_idx, seg2_idx,
- c1.pts[r1.segs[seg1_idx].begin].x,
- c1.pts[r1.segs[seg1_idx].begin].y,
- c1.pts[r1.segs[seg1_idx].end].x,
- c1.pts[r1.segs[seg1_idx].end].y,
- c2.pts[r2.segs[seg2_idx].begin].x,
- c2.pts[r2.segs[seg2_idx].begin].y,
- c2.pts[r2.segs[seg2_idx].end].x,
- c2.pts[r2.segs[seg2_idx].end].y);
- }
- p1[seg1_idx] = 1;
- std::vector<bool> &p2 = mp_problems[r2_idx];
- p2[seg2_idx] = 1;
- have_problems++;
+
+ int crosses = segs_cross(r1_idx==r2_idx,
+ c1, r1.segs[seg1_idx], c2, r2.segs[seg2_idx]);
+ if(crosses) {
+ if(VERBOSE) {
+ printf("found a crossing (still): %zd,%zd,%zd,%zd (%f,%f)-(%f,%f) (%f,%f)-(%f,%f)\n",
+ r1_idx, seg1_idx, r2_idx, seg2_idx,
+ c1.pts[r1.segs[seg1_idx].begin].x,
+ c1.pts[r1.segs[seg1_idx].begin].y,
+ c1.pts[r1.segs[seg1_idx].end].x,
+ c1.pts[r1.segs[seg1_idx].end].y,
+ c2.pts[r2.segs[seg2_idx].begin].x,
+ c2.pts[r2.segs[seg2_idx].begin].y,
+ c2.pts[r2.segs[seg2_idx].end].x,
+ c2.pts[r2.segs[seg2_idx].end].y);
}
- } // seg loop
- } // ring loop
- } // seg loop
- } // ring loop
+ p1[seg1_idx] = 1;
+ std::vector<bool> &p2 = mp_problems[r2_idx];
+ p2[seg2_idx] = 1;
+ num_problems++;
+ }
+ } // ring2/seg2 loop
+ } // seg1 loop
+ } // ring1 loop
progress += (1.0-progress)/2;
} // while problems
GDALTermProgress(1, NULL, NULL);
- if(have_problems) {
+ if(num_problems) {
printf("WARNING: Could not fix all topology problems.\n Please inspect output shapefile manually.\n");
}
}
diff --git a/src/dp.h b/src/dp.h
index a5d4ca8..feb7cf8 100644
--- a/src/dp.h
+++ b/src/dp.h
@@ -39,6 +39,13 @@ struct segment_t {
segment_t() : begin(0), end(0) { }
segment_t(size_t _begin, size_t _end) : begin(_begin), end(_end) { }
+ Bbox get_bbox(const Ring &ring) const {
+ Bbox bbox;
+ bbox.expand(ring.pts[begin]);
+ bbox.expand(ring.pts[end]);
+ return bbox;
+ }
+
size_t begin;
size_t end;
};
diff --git a/src/gdal_contrast_stretch.cc b/src/gdal_contrast_stretch.cc
index dd1a14c..90080ee 100644
--- a/src/gdal_contrast_stretch.cc
+++ b/src/gdal_contrast_stretch.cc
@@ -417,7 +417,6 @@ int main(int argc, char *argv[]) {
buf_out[band_idx].resize(block_len);
}
std::vector<uint8_t> ndv_mask(block_len);
- std::vector<uint8_t> band_mask(block_len);
for(size_t boff_y=0; boff_y<h; boff_y+=blocksize_y) {
size_t bsize_y = blocksize_y;
@@ -437,15 +436,10 @@ int main(int argc, char *argv[]) {
for(size_t band_idx=0; band_idx<dst_band_count; band_idx++) {
GDALRasterIO(src_bands[band_idx], GF_Read, boff_x, boff_y, bsize_x, bsize_y,
&buf_in[band_idx][0], bsize_x, bsize_y, GDT_Float64, 0, 0);
-
- if(band_idx == 0) {
- ndv_def.arrayCheckNdv(band_idx, &buf_in[band_idx][0], &ndv_mask[0], block_len);
- } else {
- ndv_def.arrayCheckNdv(band_idx, &buf_in[band_idx][0], &band_mask[0], block_len);
- ndv_def.aggregateMask(&ndv_mask[0], &band_mask[0], block_len);
- }
}
+ ndv_def.getNdvMask(buf_in, &ndv_mask[0], block_len);
+
for(size_t band_idx=0; band_idx<dst_band_count; band_idx++) {
double *p_in = &buf_in[band_idx][0];
uint8_t *p_out = &buf_out[band_idx][0];
@@ -517,7 +511,6 @@ std::vector<std::pair<double, double> > compute_minmax(
buf_in[band_idx].resize(block_len);
}
std::vector<uint8_t> ndv_mask(block_len);
- std::vector<uint8_t> band_mask(block_len);
std::vector<bool> got_data(band_count);
@@ -539,14 +532,10 @@ std::vector<std::pair<double, double> > compute_minmax(
for(size_t band_idx=0; band_idx<band_count; band_idx++) {
GDALRasterIO(src_bands[band_idx], GF_Read, boff_x, boff_y, bsize_x, bsize_y,
&buf_in[band_idx][0], bsize_x, bsize_y, GDT_Float64, 0, 0);
-
- if(band_idx == 0) {
- ndv_def.arrayCheckNdv(band_idx, &buf_in[band_idx][0], &ndv_mask[0], block_len);
- } else {
- ndv_def.arrayCheckNdv(band_idx, &buf_in[band_idx][0], &band_mask[0], block_len);
- ndv_def.aggregateMask(&ndv_mask[0], &band_mask[0], block_len);
- }
}
+
+ ndv_def.getNdvMask(buf_in, &ndv_mask[0], block_len);
+
for(size_t band_idx=0; band_idx<band_count; band_idx++) {
for(size_t i=0; i<block_len; i++) {
if(ndv_mask[i]) continue;
@@ -595,7 +584,6 @@ std::vector<Histogram> compute_histogram(
buf_in[band_idx].resize(block_len);
}
std::vector<uint8_t> ndv_mask(block_len);
- std::vector<uint8_t> band_mask(block_len);
std::vector<bool> first_valid_pixel(band_count, true);
@@ -617,14 +605,10 @@ std::vector<Histogram> compute_histogram(
for(size_t band_idx=0; band_idx<band_count; band_idx++) {
GDALRasterIO(src_bands[band_idx], GF_Read, boff_x, boff_y, bsize_x, bsize_y,
&buf_in[band_idx][0], bsize_x, bsize_y, GDT_Float64, 0, 0);
-
- if(band_idx == 0) {
- ndv_def.arrayCheckNdv(band_idx, &buf_in[band_idx][0], &ndv_mask[0], block_len);
- } else {
- ndv_def.arrayCheckNdv(band_idx, &buf_in[band_idx][0], &band_mask[0], block_len);
- ndv_def.aggregateMask(&ndv_mask[0], &band_mask[0], block_len);
- }
}
+
+ ndv_def.getNdvMask(buf_in, &ndv_mask[0], block_len);
+
for(size_t band_idx=0; band_idx<band_count; band_idx++) {
Histogram &hg = histograms[band_idx];
double *p = &buf_in[band_idx][0];
diff --git a/src/gdal_dem2rgb.cc b/src/gdal_dem2rgb.cc
index a6abd84..96f42a8 100644
--- a/src/gdal_dem2rgb.cc
+++ b/src/gdal_dem2rgb.cc
@@ -73,7 +73,7 @@ void usage(const std::string &cmdname) {
printf(" -alpha-overlay Generate an RGBA image that can be used as a hillshade mask\n");
printf("\n");
printf("Shading:\n");
- printf(" -exag slope_exageration Exagerate slope (default: %.1f)\n", default_slope_exageration);
+ printf(" -exag slope_exageration Exaggerate slope (default: %.1f)\n", default_slope_exageration);
printf(" -shade ambient diffuse specular_intensity specular_falloff (default: %.1f, %.1f, %.1f, %.1f)\n",
default_shade_params[0], default_shade_params[1], default_shade_params[2], default_shade_params[3]);
printf(" -lightvec sun_x sun_y sun_z (default: %.1f, %.1f, %.1f)\n",
@@ -387,9 +387,9 @@ int main(int argc, char *argv[]) {
GDALRasterIO(src_band, GF_Read, 0, 0, w, 1, &inbuf_prev[0], w, 1, GDT_Float64, 0, 0);
GDALRasterIO(src_band, GF_Read, 0, 0, w, 1, &inbuf_this[0], w, 1, GDT_Float64, 0, 0);
GDALRasterIO(src_band, GF_Read, 0, 1, w, 1, &inbuf_next[0], w, 1, GDT_Float64, 0, 0);
- ndv_def.arrayCheckNdv(0, &inbuf_prev[0], &inbuf_ndv_prev[0], w);
- ndv_def.arrayCheckNdv(0, &inbuf_this[0], &inbuf_ndv_this[0], w);
- ndv_def.arrayCheckNdv(0, &inbuf_next[0], &inbuf_ndv_next[0], w);
+ ndv_def.getNdvMask(&inbuf_prev[0], GDT_Float64, &inbuf_ndv_prev[0], w);
+ ndv_def.getNdvMask(&inbuf_this[0], GDT_Float64, &inbuf_ndv_this[0], w);
+ ndv_def.getNdvMask(&inbuf_next[0], GDT_Float64, &inbuf_ndv_next[0], w);
scale_values(&inbuf_prev[0], w, src_scale, src_offset);
scale_values(&inbuf_this[0], w, src_scale, src_offset);
scale_values(&inbuf_next[0], w, src_scale, src_offset);
@@ -408,7 +408,7 @@ int main(int argc, char *argv[]) {
int read_row = (row == h-1) ? row : row+1;
GDALRasterIO(src_band, GF_Read, 0, read_row, w, 1, &inbuf_next[0], w, 1, GDT_Float64, 0, 0);
- ndv_def.arrayCheckNdv(0, &inbuf_next[0], &inbuf_ndv_next[0], w);
+ ndv_def.getNdvMask(&inbuf_next[0], GDT_Float64, &inbuf_ndv_next[0], w);
scale_values(&inbuf_next[0], w, src_scale, src_offset);
}
if(!tex_bands.empty()) {
@@ -676,8 +676,10 @@ void compute_tierow_invaffine(
int num_cols, int row, int grid_spacing,
double *invaffine_tierow
) {
- double tiecol_left[4];
- double tiecol_right[4];
+ // this will be initialized on the first iteration, but it is set here to avoid a compiler
+ // warning.
+ double tiecol_left[4] = { 0, 0, 0, 0 };
+ double tiecol_right[4] = { 0, 0, 0, 0 };
compute_invaffine(georef, 0, row, tiecol_right);
diff --git a/src/gdal_get_projected_bounds.cc b/src/gdal_get_projected_bounds.cc
index f0512c5..03ffde7 100644
--- a/src/gdal_get_projected_bounds.cc
+++ b/src/gdal_get_projected_bounds.cc
@@ -59,7 +59,7 @@ void usage(const std::string &cmdname) {
printf(" -t_bounds_wkt <fn> File containing WKT for valid region of target SRS (optional)\n");
printf(" -s_srs <srs_def> Source SRS\n");
printf(" -t_srs <srs_def> Target SRS\n");
- printf(" -report <out.ppm> Ouput a graphical report (optional)\n");
+ printf(" -report <out.ppm> Output a graphical report (optional)\n");
printf("\nOutput is the envelope of the source region projected into the target SRS.\n");
printf("If the -t_bounds_wkt option is given it will be used as a clip mask in the\n");
printf("projected space.\n");
diff --git a/src/gdal_landsat_pansharp.cc b/src/gdal_landsat_pansharp.cc
index 9d289eb..2a721a2 100644
--- a/src/gdal_landsat_pansharp.cc
+++ b/src/gdal_landsat_pansharp.cc
@@ -32,6 +32,8 @@ This code was developed by Dan Stahlke for the Geographic Information Network of
#include <vector>
+using namespace dangdal;
+
struct ScaledBand {
ScaledBand() :
oversample(0), lo_w(0), lo_h(0), hi_w(0), hi_h(0),
diff --git a/src/gdal_merge_simple.cc b/src/gdal_merge_simple.cc
index 7dcea14..65e0736 100644
--- a/src/gdal_merge_simple.cc
+++ b/src/gdal_merge_simple.cc
@@ -30,6 +30,8 @@ This code was developed by Dan Stahlke for the Geographic Information Network of
#include "common.h"
+using namespace dangdal;
+
void copyGeoCode(GDALDatasetH dst_ds, GDALDatasetH src_ds);
void usage(const std::string &cmdname) {
diff --git a/src/gdal_merge_vrt.cc b/src/gdal_merge_vrt.cc
index 643ed59..c027639 100644
--- a/src/gdal_merge_vrt.cc
+++ b/src/gdal_merge_vrt.cc
@@ -30,6 +30,8 @@ This code was developed by Dan Stahlke for the Geographic Information Network of
#include "common.h"
+using namespace dangdal;
+
void usage(const std::string &cmdname) {
printf("Usage:\n");
printf(" %s -in <rgb.tif> -in <mask.tif> -out <out.vrt>\n", cmdname.c_str());
diff --git a/src/gdal_raw2geotiff.cc b/src/gdal_raw2geotiff.cc
index d38d694..11245e4 100644
--- a/src/gdal_raw2geotiff.cc
+++ b/src/gdal_raw2geotiff.cc
@@ -32,6 +32,8 @@ This code was developed by Dan Stahlke for the Geographic Information Network of
#include "common.h"
+using namespace dangdal;
+
void usage(const std::string &cmdname) {
printf("Usage: %s\n", cmdname.c_str());
printf("\t-wh <width> <height>\n");
diff --git a/src/gdal_trace_outline.cc b/src/gdal_trace_outline.cc
index f0556a4..a8e8694 100644
--- a/src/gdal_trace_outline.cc
+++ b/src/gdal_trace_outline.cc
@@ -26,11 +26,12 @@ This code was developed by Dan Stahlke for the Geographic Information Network of
-#include <boost/lexical_cast.hpp>
-#include <boost/foreach.hpp>
#include <map>
#include <vector>
+#include <boost/lexical_cast.hpp>
+#include <boost/foreach.hpp>
+
#include "common.h"
#include "polygon.h"
#include "polygon-rasterizer.h"
@@ -42,6 +43,7 @@ This code was developed by Dan Stahlke for the Geographic Information Network of
#include "dp.h"
#include "excursion_pincher.h"
#include "beveler.h"
+#include "raster_features.h"
#include <ogrsf_frmts.h>
#include <cpl_string.h>
@@ -151,11 +153,8 @@ struct GeomOutput {
wkt_fh(NULL),
wkb_fh(NULL),
ogr_ds(NULL),
- ogr_layer(NULL),
- class_fld_idx(-1)
- {
- for(int i=0; i<4; i++) color_fld_idx[i] = -1;
- }
+ ogr_layer(NULL)
+ { }
CoordSystem out_cs;
@@ -169,8 +168,6 @@ struct GeomOutput {
std::string ogr_fmt;
OGRDataSourceH ogr_ds;
OGRLayerH ogr_layer;
- int class_fld_idx;
- int color_fld_idx[4];
};
struct ContainingOption {
@@ -336,7 +333,6 @@ int main(int argc, char **argv) {
}
if(classify) {
- if(!ndv_def.empty()) fatal_error("-classify option is not compatible with NDV options");
if(do_invert) fatal_error("-classify option is not compatible with -invert option");
if(mask_out_fn.size()) fatal_error("-classify option is not compatible with -mask-out option");
}
@@ -347,15 +343,12 @@ int main(int argc, char **argv) {
if(!ds) fatal_error("open failed");
if(inspect_bandids.empty()) {
- size_t nbands = classify ? 1 : GDALGetRasterCount(ds);
+ size_t nbands = GDALGetRasterCount(ds);
for(size_t i=0; i<nbands; i++) inspect_bandids.push_back(i+1);
}
- // FIXME - optional NDV for classify
- if(!classify) {
- if(ndv_def.empty()) {
- ndv_def = NdvDef(ds, inspect_bandids);
- }
+ if(ndv_def.empty()) {
+ ndv_def = NdvDef(ds, inspect_bandids);
}
CPLPushErrorHandler(CPLQuietErrorHandler);
@@ -378,29 +371,17 @@ int main(int argc, char **argv) {
do_pinch_excursions ? PLOT_PINCH : PLOT_CONTOURS);
}
- std::vector<uint8_t> raster;
- BitGrid mask(0, 0);
- uint8_t usage_array[256];
- GDALColorTableH color_table = NULL;
- if(classify) {
- if(inspect_bandids.size() != 1) {
- fatal_error("only one band may be used in classify mode");
- }
-
- std::vector<uint8_t> tmp_raster = read_dataset_8bit(ds, inspect_bandids[0], usage_array, dbuf);
- std::swap(raster, tmp_raster);
+ // only used if classify==true, but doesn't hurt otherwise
+ const FeatureInterpreter feature_interp(ds, inspect_bandids);
- GDALRasterBandH band = GDALGetRasterBand(ds, inspect_bandids[0]);
- if(GDALGetRasterColorInterpretation(band) == GCI_PaletteIndex) {
- color_table = GDALGetRasterColorTable(band);
- }
- } else {
- mask = get_bitgrid_for_dataset(ds, inspect_bandids, ndv_def, dbuf);
+ FeatureBitmap *features_bitmap = NULL;
+ if(classify) {
+ features_bitmap = FeatureBitmap::from_raster(ds, inspect_bandids, ndv_def, dbuf);
}
for(size_t go_idx=0; go_idx<geom_outputs.size(); go_idx++) {
GeomOutput &go = geom_outputs[go_idx];
-
+
if(go.wkt_fn.size()) {
go.wkt_fh = fopen(go.wkt_fn.c_str(), "w");
if(!go.wkt_fh) fatal_error("cannot open output file for WKT");
@@ -416,7 +397,8 @@ int main(int argc, char **argv) {
OGRSFDriverH ogr_driver = OGRGetDriverByName(go.ogr_fmt.c_str());
if(!ogr_driver) fatal_error("cannot get OGR driver (%s)", go.ogr_fmt.c_str());
go.ogr_ds = OGR_Dr_CreateDataSource(ogr_driver, go.ogr_fn.c_str(), NULL);
- if(!go.ogr_ds) fatal_error("cannot create OGR data source");
+ if(!go.ogr_ds) fatal_error(
+ "cannot create OGR data source (does output file already exist?)");
std::string layer_name = go.ogr_fn;
@@ -432,50 +414,42 @@ int main(int argc, char **argv) {
if(!go.ogr_layer) fatal_error("cannot create OGR layer");
if(classify) {
- OGRFieldDefnH fld = OGR_Fld_Create("value", OFTInteger);
- OGR_Fld_SetWidth(fld, 4);
- OGR_L_CreateField(go.ogr_layer, fld, TRUE);
- go.class_fld_idx = OGR_FD_GetFieldIndex(OGR_L_GetLayerDefn(go.ogr_layer), "value");
-
- if(color_table) {
- const char *names[4] = { "c1", "c2", "c3", "c4" };
- for(int i=0; i<4; i++) {
- fld = OGR_Fld_Create(names[i], OFTInteger);
- OGR_Fld_SetWidth(fld, 4);
- OGR_L_CreateField(go.ogr_layer, fld, TRUE);
- go.color_fld_idx[i] = OGR_FD_GetFieldIndex(OGR_L_GetLayerDefn(go.ogr_layer), names[i]);
- }
- }
+ feature_interp.create_ogr_fields(go.ogr_layer);
}
}
}
int num_shapes_written = 0;
- for(int class_id=0; class_id<256; class_id++) {
- const GDALColorEntry *color = NULL;
- if(classify) {
- if(!usage_array[class_id]) continue;
- printf("\nFeature class %d\n", class_id);
+ std::map<FeatureRawVal, FeatureBitmap::Index> features_list;
+ if(classify) {
+ features_list = features_bitmap->feature_table();
+ } else {
+ // If we are not running in feature classify mode, then just add an arbitrary element
+ // to this list so that the loop below runs once.
+ features_list[FeatureRawVal()];
+ }
- if(color_table) {
- color = GDALGetColorEntry(color_table, class_id);
- if(color) printf(" Color=%d,%d,%d,%d\n",
- color->c1, color->c2, color->c3, color->c4);
+ typedef std::map<FeatureRawVal, FeatureBitmap::Index>::value_type feature_pair_t;
+ size_t feature_idx = 0;
+ BOOST_FOREACH(const feature_pair_t &feature, features_list) {
+ Mpoly feature_poly;
+ {
+ BitGrid mask(0, 0);
+ if(classify) {
+ printf("\nTracing feature %s (%zd of %zd)\n",
+ feature_interp.pixel_to_string(feature.first).c_str(),
+ (++feature_idx), features_list.size());
+ mask = features_bitmap->get_mask_for_feature(feature.second);
+ } else {
+ printf("Reading raster.\n");
+ mask = get_bitgrid_for_dataset(ds, inspect_bandids, ndv_def, dbuf);
}
- mask = get_bitgrid_for_8bit_raster(georef.w, georef.h,
- &raster[0], (uint8_t)class_id);
- } else {
- if(class_id != 0) continue;
- }
-
- if(do_invert) {
- mask.invert();
- }
+ if(do_invert) mask.invert();
+ if(do_erosion) mask.erode();
- if(do_erosion) {
- mask.erode();
+ feature_poly = trace_mask(mask, georef.w, georef.h, min_ring_area, trace_no_donuts);
}
if(!containing_options.empty()) {
@@ -499,9 +473,6 @@ int main(int argc, char **argv) {
trace_no_donuts = 1;
}
- Mpoly feature_poly = trace_mask(mask, georef.w, georef.h, min_ring_area, trace_no_donuts);
- mask = BitGrid(0, 0); // free some memory
-
if(VERBOSE) {
size_t num_inner = 0, num_outer = 0, total_pts = 0;
for(size_t r_idx=0; r_idx<feature_poly.rings.size(); r_idx++) {
@@ -609,17 +580,11 @@ int main(int argc, char **argv) {
if(go.ogr_ds) {
OGRFeatureH ogr_feat = OGR_F_Create(OGR_L_GetLayerDefn(go.ogr_layer));
- if(go.class_fld_idx >= 0) OGR_F_SetFieldInteger(ogr_feat, go.class_fld_idx, class_id);
- if(color) {
- if(go.color_fld_idx[0] >= 0) OGR_F_SetFieldInteger(
- ogr_feat, go.color_fld_idx[0], color->c1);
- if(go.color_fld_idx[1] >= 0) OGR_F_SetFieldInteger(
- ogr_feat, go.color_fld_idx[1], color->c2);
- if(go.color_fld_idx[2] >= 0) OGR_F_SetFieldInteger(
- ogr_feat, go.color_fld_idx[2], color->c3);
- if(go.color_fld_idx[3] >= 0) OGR_F_SetFieldInteger(
- ogr_feat, go.color_fld_idx[3], color->c4);
+
+ if(classify) {
+ feature_interp.set_ogr_fields(go.ogr_layer, ogr_feat, feature.first);
}
+
OGR_F_SetGeometryDirectly(ogr_feat, ogr_geom); // assumes ownership of geom
OGR_L_CreateFeature(go.ogr_layer, ogr_feat);
OGR_F_Destroy(ogr_feat);
@@ -636,6 +601,8 @@ int main(int argc, char **argv) {
printf("\n");
+ delete(features_bitmap);
+
for(size_t go_idx=0; go_idx<geom_outputs.size(); go_idx++) {
GeomOutput &go = geom_outputs[go_idx];
if(go.wkt_fh) fclose(go.wkt_fh);
diff --git a/src/mask-tracer.cc b/src/mask-tracer.cc
index dfa79c7..6532c78 100644
--- a/src/mask-tracer.cc
+++ b/src/mask-tracer.cc
@@ -104,10 +104,10 @@ static inline pixquad_t get_quad(const BitGrid &mask, int x, int y, bool select_
// 1 2
// 8 4
pixquad_t quad =
- (mask(x-1, y-1) ? 1 : 0) +
- (mask(x , y-1) ? 2 : 0) +
- (mask(x , y ) ? 4 : 0) +
- (mask(x-1, y ) ? 8 : 0);
+ (mask.get(x-1, y-1, 0) ? 1 : 0) +
+ (mask.get(x , y-1, 0) ? 2 : 0) +
+ (mask.get(x , y , 0) ? 4 : 0) +
+ (mask.get(x-1, y , 0) ? 8 : 0);
if(!select_color) quad ^= 0xf;
return quad;
}
@@ -166,31 +166,20 @@ int initial_x, int initial_y, bool select_color) {
}
static int recursive_trace(BitGrid &mask, size_t w, size_t h,
-const Ring &bounds, int depth, Mpoly &out_poly, int parent_id,
+const Ring &bounding_ring, int depth, Mpoly &out_poly, int parent_id,
int64_t min_area, bool no_donuts) {
//printf("recursive_trace enter: depth=%d\n", depth);
bool select_color = !(depth & 1);
- int bound_left, bound_right;
- int bound_top, bound_bottom;
- bound_left = bound_right = (int)bounds.pts[0].x;
- bound_top = bound_bottom = (int)bounds.pts[0].y;
- for(size_t v_idx=0; v_idx<bounds.pts.size(); v_idx++) {
- int x = (int)bounds.pts[v_idx].x;
- int y = (int)bounds.pts[v_idx].y;
- if(x < bound_left) bound_left = x;
- if(x > bound_right) bound_right = x;
- if(y < bound_top) bound_top = y;
- if(y > bound_bottom) bound_bottom = y;
- }
+ Bbox bounding_bbox = bounding_ring.getBbox();
Mpoly bounds_mp;
- bounds_mp.rings.push_back(bounds);
+ bounds_mp.rings.push_back(bounding_ring);
- std::vector<row_crossings_t> crossings =
- get_row_crossings(bounds_mp, bound_top, bound_bottom-bound_top);
- assert(crossings.size() == size_t(bound_bottom-bound_top));
+ std::vector<row_crossings_t> crossings =
+ get_row_crossings(bounds_mp, bounding_bbox.min_y, bounding_bbox.height());
+ assert(crossings.size() == size_t(bounding_bbox.height()));
int skip_this = min_area && (compute_area(crossings) < min_area);
int skip_child = skip_this || (depth && no_donuts);
@@ -200,63 +189,48 @@ int64_t min_area, bool no_donuts) {
}
if(!skip_child) {
- for(int y=bound_top+1; y<bound_bottom; y++) {
+ for(int y=bounding_bbox.min_y+1; y<bounding_bbox.max_y; y++) {
if(!depth) {
- GDALTermProgress((double)y/(double)(bound_bottom-bound_top-1), NULL, NULL);
+ GDALTermProgress((double)y/(double)(bounding_bbox.height()-1), NULL, NULL);
}
// make sure the range (y-1,y)*(x-1,x) is in bounds
row_crossings_t cross_both = crossings_intersection(
- crossings[y-bound_top-1], crossings[y-bound_top]);
+ crossings[y-bounding_bbox.min_y-1], crossings[y-bounding_bbox.min_y]);
for(size_t cidx=0; cidx<cross_both.size()/2; cidx++) {
// make sure the range (y-1,y)*(x-1,x) is in bounds
int from = 1+cross_both[cidx*2 ];
int to = cross_both[cidx*2+1];
for(int x=from; x<to; x++) {
- //pixquad_t quad = get_quad(mask, x, y, select_color);
- //int is_seed = (quad != 0 && quad != 0xf);
- int is_seed;
- if(select_color) {
- is_seed =
- mask(x-1, y-1) ||
- mask(x , y-1) ||
- mask(x-1, y ) ||
- mask(x , y );
- } else {
- is_seed =
- (!mask(x-1, y-1)) ||
- (!mask(x , y-1)) ||
- (!mask(x-1, y )) ||
- (!mask(x , y ));
- }
+ pixquad_t quad = get_quad(mask, x, y, select_color);
+ int is_seed = (quad != 0);
if(is_seed) {
Ring r = trace_single_mpoly(mask, w, h, x, y, select_color);
r.parent_id = parent_id;
r.is_hole = depth % 2;
- //r.parent_id = -1;
- //r.is_hole = 0;
-
size_t outer_ring_id = out_poly.rings.size();
out_poly.rings.push_back(r);
int was_skip = recursive_trace(
mask, w, h, r, depth+1, out_poly, outer_ring_id,
min_area, no_donuts);
+
if(was_skip) {
out_poly.rings.pop_back();
}
}
- }
+ }
}
}
}
- for(int y=bound_top; y<bound_bottom; y++) {
- const row_crossings_t &r = crossings[y-bound_top];
- if(depth>0) {
+ if(depth>0) {
+ // erase this polygon from the raster by filling it with select_color
+ for(int y=bounding_bbox.min_y; y<bounding_bbox.max_y; y++) {
+ const row_crossings_t &r = crossings[y-bounding_bbox.min_y];
for(size_t cidx=0; cidx<r.size()/2; cidx++) {
int from = r[cidx*2 ];
int to = r[cidx*2+1];
diff --git a/src/mask.cc b/src/mask.cc
index 268bd38..ff28d65 100644
--- a/src/mask.cc
+++ b/src/mask.cc
@@ -27,230 +27,145 @@ This code was developed by Dan Stahlke for the Geographic Information Network of
#include <vector>
+#include <algorithm>
+
+#include <boost/foreach.hpp>
#include "common.h"
#include "mask.h"
#include "polygon.h"
#include "debugplot.h"
#include "ndv.h"
+#include "datatype_conversion.h"
namespace dangdal {
-std::vector<uint8_t> read_dataset_8bit(GDALDatasetH ds, int band_idx, uint8_t *usage_array, DebugPlot *dbuf) {
- for(int i=0; i<256; i++) usage_array[i] = 0;
+BitGrid get_bitgrid_for_dataset(
+ GDALDatasetH ds, const std::vector<size_t> &band_ids,
+ const NdvDef &ndv_def, DebugPlot *dbuf
+) {
+ assert(!band_ids.empty());
size_t w = GDALGetRasterXSize(ds);
size_t h = GDALGetRasterYSize(ds);
- int band_count = GDALGetRasterCount(ds);
- if(VERBOSE) printf("input is %zd x %zd x %d\n", w, h, band_count);
-
- if(band_idx < 1 || band_idx > band_count) fatal_error("bandid out of range");
+ size_t band_count = GDALGetRasterCount(ds);
+ if(VERBOSE) printf("input is %zd x %zd x %zd\n", w, h, band_count);
- GDALRasterBandH band = GDALGetRasterBand(ds, band_idx);
+ std::vector<GDALRasterBandH> bands;
+ std::vector<GDALDataType> datatypes;
+ BOOST_FOREACH(const size_t band_id, band_ids) {
+ if(VERBOSE) printf("opening band %zd\n", band_id);
+ GDALRasterBandH band = GDALGetRasterBand(ds, band_id);
+ if(!band) fatal_error("Could not open band %zd.", band_id);
+ bands.push_back(band);
+ GDALDataType dt = GDALGetRasterDataType(band);
+ datatypes.push_back(dt);
+ }
int blocksize_x_int, blocksize_y_int;
- GDALGetBlockSize(band, &blocksize_x_int, &blocksize_y_int);
+ GDALGetBlockSize(bands[0], &blocksize_x_int, &blocksize_y_int);
+ // Out of laziness, I am hoping that images always have the same block size for each band.
+ BOOST_FOREACH(const GDALRasterBandH band, bands) {
+ int bx, by;
+ GDALGetBlockSize(band, &bx, &by);
+ if(bx != blocksize_x_int || by != blocksize_y_int) {
+ fatal_error(
+ "Bands have different block sizes. Not currently implemented. Please contact developer");
+ }
+ }
size_t blocksize_x = blocksize_x_int;
size_t blocksize_y = blocksize_y_int;
+ size_t blocksize_xy = blocksize_x * blocksize_y;
- GDALDataType gdt = GDALGetRasterDataType(band);
- if(gdt != GDT_Byte) {
- printf("Warning: input is not of type Byte, there may be loss while downsampling!\n");
+ std::vector<std::vector<uint8_t> > band_buf(bands.size());
+ for(size_t i=0; i<bands.size(); i++) {
+ size_t dt_size = GDALGetDataTypeSize(datatypes[i]) / 8;
+ size_t num_bytes = blocksize_xy * dt_size;
+ band_buf[i].resize(num_bytes);
}
- if(VERBOSE) printf("band %d: block size = %zd,%zd\n",
- band_idx, blocksize_x, blocksize_y);
+ std::vector<uint8_t> block_mask(blocksize_xy);
+
+ BitGrid mask(w, h);
+ mask.zero();
+
+ size_t num_valid = 0;
+ size_t num_ndv = 0;
- printf("Reading one band of size %zd x %zd\n", w, h);
+ printf("Reading input...\n");
- std::vector<uint8_t> outbuf(w*h);
- std::vector<uint8_t> inbuf(blocksize_x*blocksize_y);
- for(size_t boff_y=0; boff_y<h; boff_y+=blocksize_y) {
+ size_t num_blocks_x = (w + blocksize_x - 1) / blocksize_x;
+ size_t num_blocks_y = (h + blocksize_y - 1) / blocksize_y;
+ for(size_t block_y=0; block_y<num_blocks_y; block_y++) {
+ size_t boff_y = blocksize_y * block_y;
size_t bsize_y = blocksize_y;
if(bsize_y + boff_y > h) bsize_y = h - boff_y;
- for(size_t boff_x=0; boff_x<w; boff_x+=blocksize_x) {
+ for(size_t block_x=0; block_x<num_blocks_x; block_x++) {
+ size_t boff_x = blocksize_x * block_x;
size_t bsize_x = blocksize_x;
if(bsize_x + boff_x > w) bsize_x = w - boff_x;
- double progress =
+ double progress =
double(
boff_y * w +
boff_x * bsize_y
) / (w * h);
GDALTermProgress(progress, NULL, NULL);
- GDALRasterIO(band, GF_Read, boff_x, boff_y, bsize_x, bsize_y,
- &inbuf[0], bsize_x, bsize_y, GDT_Byte, 0, 0);
-
- uint8_t *p_in = &inbuf[0];
- for(size_t j=0; j<bsize_y; j++) {
- size_t y = j + boff_y;
- bool is_dbuf_stride_y = dbuf && ((y % dbuf->stride_y) == 0);
- uint8_t *p_out = &outbuf[w*y + boff_x];
- for(size_t i=0; i<bsize_x; i++) {
- uint8_t val = *(p_in++);
- *(p_out++) = val;
- usage_array[val] = 1;
-
- bool is_dbuf_stride = is_dbuf_stride_y && ((i % dbuf->stride_x) == 0);
- if(is_dbuf_stride) {
- size_t x = i + boff_x;
- uint8_t db_v = 50 + val/3;
- if(db_v < 50) db_v = 50;
- if(db_v > 254) db_v = 254;
- uint8_t r = (uint8_t)(db_v*.75);
- dbuf->plotPoint(x, y, r, db_v, db_v);
- }
- }
+ for(size_t band_idx=0; band_idx<bands.size(); band_idx++) {
+ GDALReadBlock(bands[band_idx], block_x, block_y, &band_buf[band_idx][0]);
}
- }
- }
-
- GDALTermProgress(1, NULL, NULL);
-
- return outbuf;
-}
-
-BitGrid get_bitgrid_for_dataset(
- GDALDatasetH ds, const std::vector<size_t> &bandlist,
- const NdvDef &ndv_def, DebugPlot *dbuf
-) {
- size_t w = GDALGetRasterXSize(ds);
- size_t h = GDALGetRasterYSize(ds);
- size_t band_count = GDALGetRasterCount(ds);
- if(VERBOSE) printf("input is %zd x %zd x %zd\n", w, h, band_count);
-
- BitGrid mask(w, h);
- mask.zero();
-
- printf("Reading %zd bands of size %zd x %zd\n", bandlist.size(), w, h);
-
- for(size_t bandlist_idx=0; bandlist_idx<bandlist.size(); bandlist_idx++) {
- size_t band_idx = bandlist[bandlist_idx];
- if(band_idx < 1 || band_idx > band_count) fatal_error("bandid out of range");
-
- GDALRasterBandH band = GDALGetRasterBand(ds, band_idx);
-
- int blocksize_x_int, blocksize_y_int;
- GDALGetBlockSize(band, &blocksize_x_int, &blocksize_y_int);
- size_t blocksize_x = blocksize_x_int;
- size_t blocksize_y = blocksize_y_int;
-
- // gain speed in common 8-bit case
- GDALDataType gdt = GDALGetRasterDataType(band);
- bool use_8bit = (gdt == GDT_Byte);
-
- if(VERBOSE) printf("band %zd: block size = %zd,%zd, use_8bit=%d\n",
- band_idx, blocksize_x, blocksize_y, use_8bit?1:0);
-
- std::vector<uint8_t> block_buf_8bit;
- std::vector<double> block_buf_dbl;
- if(use_8bit) {
- block_buf_8bit.resize(blocksize_x*blocksize_y);
- } else {
- block_buf_dbl.resize(blocksize_x*blocksize_y);
- }
-
- std::vector<uint8_t> row_ndv(blocksize_x);
-
- for(size_t boff_y=0; boff_y<h; boff_y+=blocksize_y) {
- size_t bsize_y = blocksize_y;
- if(bsize_y + boff_y > h) bsize_y = h - boff_y;
- for(size_t boff_x=0; boff_x<w; boff_x+=blocksize_x) {
- size_t bsize_x = blocksize_x;
- if(bsize_x + boff_x > w) bsize_x = w - boff_x;
-
- double progress =
- double(
- bandlist_idx * w * h +
- boff_y * w +
- boff_x * bsize_y
- ) / (bandlist.size() * w * h);
- GDALTermProgress(progress, NULL, NULL);
-
- if(use_8bit) {
- GDALRasterIO(band, GF_Read, boff_x, boff_y, bsize_x, bsize_y,
- &block_buf_8bit[0], bsize_x, bsize_y, GDT_Byte, 0, 0);
- } else {
- GDALRasterIO(band, GF_Read, boff_x, boff_y, bsize_x, bsize_y,
- &block_buf_dbl[0], bsize_x, bsize_y, GDT_Float64, 0, 0);
- }
-
- double *p_dbl = NULL;
- uint8_t *p_8bit = NULL;
- if(use_8bit) {
- p_8bit = &block_buf_8bit[0];
- } else {
- p_dbl = &block_buf_dbl[0];
- }
- for(size_t j=0; j<bsize_y; j++) {
- size_t y = j + boff_y;
- bool is_dbuf_stride_y = dbuf && (bandlist_idx==0) && ((y % dbuf->stride_y) == 0);
+ ndv_def.getNdvMask(band_buf, datatypes, &block_mask[0], blocksize_xy);
- if(p_dbl) {
- ndv_def.arrayCheckNdv(bandlist_idx, p_dbl, &row_ndv[0], bsize_x);
+ for(size_t sub_y=0; sub_y<bsize_y; sub_y++) {
+ size_t y = sub_y + boff_y;
+ bool is_dbuf_stride_y = dbuf && ((y % dbuf->stride_y) == 0);
+ for(size_t sub_x=0; sub_x<bsize_x; sub_x++) {
+ size_t x = sub_x + boff_x;
+ bool is_dbuf_stride = is_dbuf_stride_y && ((sub_x % dbuf->stride_x) == 0);
+
+ bool is_ndv = block_mask[sub_y*blocksize_x + sub_x];
+ mask.set(x, y, !is_ndv);
+ if(is_ndv) {
+ num_ndv++;
} else {
- ndv_def.arrayCheckNdv(bandlist_idx, p_8bit, &row_ndv[0], bsize_x);
+ num_valid++;
}
- for(size_t i=0; i<bsize_x; i++) {
- bool is_dbuf_stride = is_dbuf_stride_y && ((i % dbuf->stride_x) == 0);
- if(is_dbuf_stride) {
- int val = use_8bit ? (int)(*p_8bit) : (int)(*p_dbl);
- size_t x = i + boff_x;
- int db_v = 50 + val/3;
- if(db_v < 50) db_v = 50;
- if(db_v > 254) db_v = 254;
- uint8_t r = (uint8_t)(db_v*.75);
- dbuf->plotPoint(x, y, r, (uint8_t)db_v, (uint8_t)db_v);
- }
-
- if(use_8bit) p_8bit++;
- else p_dbl++;
- }
-
- if(!bandlist_idx) {
- for(size_t i=0; i<bsize_x; i++) {
- mask.set(boff_x+i, y, !row_ndv[i]);
- }
- } else if(ndv_def.isInvert()) {
- for(size_t i=0; i<bsize_x; i++) {
- if(row_ndv[i]) mask.set(boff_x+i, y, false);
- }
- } else {
- for(size_t i=0; i<bsize_x; i++) {
- if(!row_ndv[i]) mask.set(boff_x+i, y, true);
+ if(is_dbuf_stride) {
+ uint8_t val[3] = { 0, 0, 0 };
+ if(!is_ndv) {
+ for(size_t rgb_idx=0; rgb_idx<3; rgb_idx++) {
+ size_t band_idx = std::min(rgb_idx, bands.size()-1);
+ double dbl_val = gdal_scalar_to_double(
+ &band_buf[band_idx][sub_y*blocksize_x + sub_x], datatypes[band_idx]);
+ // valid pixels have texture of the image, but with a cyanish hue
+ if(rgb_idx==0) {
+ val[rgb_idx] = std::max(0.0, std::min(127.0, dbl_val*0.5));
+ } else {
+ val[rgb_idx] = std::max(64.0, std::min(191.0, dbl_val*0.5+64));
+ }
+ }
}
+ dbuf->plotPoint(x, y, val[0], val[1], val[2]);
+
+ // Old color scheme:
+ //int val = gdal_scalar_to_int32(
+ // &band_buf[0][sub_y*blocksize_x + sub_x], datatypes[0]);
+ //int db_v = 50 + val/3;
+ //if(db_v < 50) db_v = 50;
+ //if(db_v > 254) db_v = 254;
+ //uint8_t r = (uint8_t)(db_v*.75);
+ //dbuf->plotPoint(x, y, r, (uint8_t)db_v, (uint8_t)db_v);
}
}
}
}
}
- if(dbuf) {
- for(size_t y=0; y<h; y+=dbuf->stride_y) {
- for(size_t x=0; x<w; x+=dbuf->stride_x) {
- if(!mask(x, y)) {
- dbuf->plotPoint(x, y, 0, 0, 0);
- }
- }
- }
- }
-
GDALTermProgress(1, NULL, NULL);
- return mask;
-}
-
-BitGrid get_bitgrid_for_8bit_raster(size_t w, size_t h, const uint8_t *raster, uint8_t wanted) {
- BitGrid mask(w, h);
-
- const uint8_t *p = raster;
- for(size_t y=0; y<h; y++) {
- for(size_t x=0; x<w; x++) {
- mask.set(x, y, *(p++) == wanted);
- }
- }
+ printf("Found %zd valid and %zd NDV pixels.\n", num_valid, num_ndv);
return mask;
}
@@ -261,14 +176,14 @@ void BitGrid::erode() { // FIXME - untested
bool *rowl = new bool[w];
for(int i=0; i<w; i++) {
rowm[i] = 0;
- rowl[i] = get(i, 0);
+ rowl[i] = get(i, 0, 0);
}
for(int y=0; y<h; y++) {
bool *tmp = rowu;
rowu = rowm; rowm = rowl; rowl = tmp;
for(int i=0; i<w; i++) {
- rowl[i] = get(i, y+1);
+ rowl[i] = get(i, y+1, 0);
}
bool ul = 0, um = rowu[0];
diff --git a/src/mask.h b/src/mask.h
index 28e1c78..3ec6d18 100644
--- a/src/mask.h
+++ b/src/mask.h
@@ -40,66 +40,80 @@ This code was developed by Dan Stahlke for the Geographic Information Network of
namespace dangdal {
-class BitGrid {
+template <typename T>
+class GridArray {
public:
- BitGrid(int _w, int _h) :
+ GridArray(int _w, int _h) :
w(_w), h(_h),
- arrlen((size_t(w)*h+7)/8),
- grid(arrlen)
+ grid(w*h)
{ }
// default dtor, copy, assign are OK
public:
- inline bool operator()(int x, int y) const { return get(x, y); }
+ typename std::vector<T>::const_reference operator()(int x, int y) const {
+ // out-of-bounds used to be okay and return 'false' but not now
+ // FIXME - make sure that is okay
+ assert(x>=0 && y>=0 && x<w && y<h);
+ return grid[size_t(y)*w + x];
+ }
+
+ typename std::vector<T>::reference operator()(int x, int y) {
+ // out-of-bounds used to be okay and return 'false' but not now
+ // FIXME - make sure that is okay
+ assert(x>=0 && y>=0 && x<w && y<h);
+ return grid[size_t(y)*w + x];
+ }
- inline bool get(int x, int y) const {
- // out-of-bounds is OK and returns false
+ T get(
+ int x, int y, const T &default_val
+ ) const {
if(x>=0 && y>=0 && x<w && y<h) {
- size_t p = size_t(y)*w + x;
- return grid[p/8] & (1 << (p&7));
+ return (*this)(x, y);
} else {
- return false;
+ return default_val;
}
}
- inline void set(int x, int y, bool val) {
- assert(x>=0 && y>=0 && x<w && y<h);
+ // FIXME - deprecate
+ const typename std::vector<T>::const_reference get(int x, int y) const {
+ return (*this)(x, y);
+ }
- size_t p = size_t(y)*w + x;
- if(val) {
- grid[p/8] |= (1 << (p&7));
- } else {
- grid[p/8] &= ~(1 << (p&7));
- }
+ // FIXME - deprecate
+ void set(int x, int y, const T &val) {
+ (*this)(x, y) = val;
}
void zero() {
- for(size_t i=0; i<arrlen; i++) {
+ for(size_t i=0; i<grid.size(); i++) {
grid[i] = 0;
}
}
+protected:
+ int w, h;
+ std::vector<T> grid;
+};
+
+class BitGrid : public GridArray<bool> {
+public:
+ BitGrid(int _w, int _h) : GridArray<bool>(_w, _h) { }
+
void invert() {
- for(size_t i=0; i<arrlen; i++) {
- grid[i] = ~grid[i];
+ for(size_t i=0; i<grid.size(); i++) {
+ grid[i] = !grid[i];
}
}
void erode();
Vertex centroid();
-
-private:
- int w, h;
- size_t arrlen;
- std::vector<uint8_t> grid;
};
-std::vector<uint8_t> read_dataset_8bit(GDALDatasetH ds, int band_idx, uint8_t *usage_array, DebugPlot *dbuf);
-BitGrid get_bitgrid_for_dataset(GDALDatasetH ds, const std::vector<size_t> &bandlist,
+// Returns a BitGrid with 'true' values correspond to valid (not ndv) pixels.
+BitGrid get_bitgrid_for_dataset(GDALDatasetH ds, const std::vector<size_t> &bandlist,
const NdvDef &ndv_def, DebugPlot *dbuf);
-BitGrid get_bitgrid_for_8bit_raster(size_t w, size_t h, const uint8_t *raster, uint8_t wanted);
} // namespace dangdal
diff --git a/src/ndv.cc b/src/ndv.cc
index 3f9e901..43507b1 100644
--- a/src/ndv.cc
+++ b/src/ndv.cc
@@ -27,12 +27,16 @@ This code was developed by Dan Stahlke for the Geographic Information Network of
#include <algorithm>
+#include <cstring>
+#include <cassert>
+#include <boost/foreach.hpp>
#include <boost/lexical_cast.hpp>
#include <boost/tokenizer.hpp>
#include "common.h"
#include "ndv.h"
+#include "datatype_conversion.h"
void usage(const std::string &cmdname); // externally defined
@@ -44,14 +48,19 @@ void NdvDef::printUsage() {
" -ndv val Set a no-data value\n"
" -ndv 'val val ...' Set a no-data value using all input bands\n"
" -ndv 'min..max min..max ...' Set a range of no-data values\n"
-" (-Inf and Inf are allowed)\n"
+" (-Inf and Inf are allowed; '*' == '-Inf..Inf')\n"
" -valid-range 'min..max min..max ...' Set a range of valid data values\n"
);
}
-NdvInterval::NdvInterval(const std::string &s) {
+NdvInterval::NdvInterval(const std::string &s_in) {
+ // copy the string in case we want to change it
+ std::string s = s_in;
+
if(VERBOSE >= 2) printf("minmax [%s]\n", s.c_str());
+ if(s == "*") s = "-Inf..Inf";
+
double min, max;
size_t delim = s.find("..");
try {
@@ -60,7 +69,7 @@ NdvInterval::NdvInterval(const std::string &s) {
} else {
std::string s1 = s.substr(0, delim);
std::string s2 = s.substr(delim+2);
- // FIXME - allow -Inf, Inf
+ // This is capable of interpreting -Inf and Inf
min = boost::lexical_cast<double>(s1);
max = boost::lexical_cast<double>(s2);
}
@@ -167,107 +176,76 @@ void NdvDef::debugPrint() const {
printf("=== end NDV\n");
}
-template<class T>
-void flagMatches(
- const NdvInterval &range,
- const T *in_data,
- uint8_t *mask_out,
- size_t nsamps
-) {
- for(size_t i=0; i<nsamps; i++) {
- T val = in_data[i];
- if(range.contains(val)) mask_out[i] = 1;
- }
+template <typename T>
+static inline bool contains_templated(const NdvInterval *interval, const void *p) {
+ T val = *(reinterpret_cast<const T *>(p));
+ return interval->contains(val);
}
-template<>
-void flagMatches<uint8_t>(
- const NdvInterval &range,
- const uint8_t *in_data,
- uint8_t *mask_out,
- size_t nsamps
-) {
- uint8_t min_byte = (uint8_t)std::max(ceil (range.first ), 0.0);
- uint8_t max_byte = (uint8_t)std::min(floor(range.second), 255.0);
- for(size_t i=0; i<nsamps; i++) {
- uint8_t v = in_data[i];
- uint8_t match = (v >= min_byte) && (v <= max_byte);
- if(match) mask_out[i] = 1;
- }
+inline bool NdvInterval::contains(const void *p, GDALDataType dt) const {
+ return DANGDAL_RUNTIME_TEMPLATE(dt, contains_templated, this, p);
}
-template<class T>
-void flagNaN(
- const T *in_data,
- uint8_t *mask_out,
- size_t nsamps
-) {
- for(size_t i=0; i<nsamps; i++) {
- if(std::isnan(in_data[i])) mask_out[i] = 1;
+void NdvDef::getNdvMask(
+ const std::vector<const void *> &bands,
+ const std::vector<GDALDataType> &dt_list,
+ uint8_t *mask_out, size_t num_pixels
+) const {
+ assert(bands.size() == dt_list.size());
+ // use uint8_t to allow incrementing the pointer by bytes
+ std::vector<const uint8_t *> in_p;
+ std::vector<size_t> dt_sizes;
+ for(size_t i=0; i<bands.size(); i++) {
+ in_p.push_back(reinterpret_cast<const uint8_t *>(bands[i]));
+ dt_sizes.push_back(GDALGetDataTypeSize(dt_list[i]) / 8);
}
-}
-template<>
-void flagNaN<uint8_t>(
- const uint8_t *in_data __attribute__((unused)),
- uint8_t *mask_out __attribute__((unused)),
- size_t nsamps __attribute__((unused))
-) { } // no-op
-
-template<class T>
-void NdvDef::arrayCheckNdv(
- size_t band, const T *in_data,
- uint8_t *mask_out, size_t nsamps
-) const {
- for(size_t i=0; i<nsamps; i++) mask_out[i] = 0;
- for(size_t slab_idx=0; slab_idx<slabs.size(); slab_idx++) {
- const NdvSlab &slab = slabs[slab_idx];
- NdvInterval range;
- if(band > 0 && slab.range_by_band.size() == 1) {
- // if only a single range is defined, use it for all bands
- range = slab.range_by_band[0];
- } else if(band < slab.range_by_band.size()) {
- range = slab.range_by_band[band];
- } else {
- fatal_error("wrong number of bands in NDV def");
- }
- flagMatches(range, in_data, mask_out, nsamps);
+ BOOST_FOREACH(const NdvSlab &slab, slabs) {
+ size_t num_intervals = slab.range_by_band.size();
+ assert(num_intervals == bands.size() || num_intervals == 1);
}
- if(invert) {
- for(size_t i=0; i<nsamps; i++) {
- mask_out[i] = mask_out[i] ? 0 : 1;
+
+ memset(mask_out, 0, num_pixels);
+ for(size_t pix_idx=0; pix_idx<num_pixels; pix_idx++) {
+ mask_out[pix_idx] = 0;
+
+ // A NaN value on any band makes this pixel NDV.
+ for(size_t i=0; i<bands.size(); i++) {
+ if(gdal_scalar_pointer_isnan(in_p[i], dt_list[i])) {
+ mask_out[pix_idx] = 1;
+ }
}
- }
- //printf("XXX %d\n", mask_out[0]?1:0);
- flagNaN(in_data, mask_out, nsamps);
-}
-void NdvDef::aggregateMask(
- uint8_t *total_mask,
- const uint8_t *band_mask,
- size_t nsamps
-) const {
- if(invert) {
- // pixel is valid only if all bands are within valid range
- for(size_t i=0; i<nsamps; i++) {
- if(band_mask[i]) total_mask[i] = 1;
+ BOOST_FOREACH(const NdvSlab &slab, slabs) {
+ bool all_match = true;
+ for(size_t i=0; i<bands.size(); i++) {
+ size_t num_intervals = slab.range_by_band.size();
+ // if only one interval is given, use it for all bands
+ size_t j = num_intervals==1 ? 0 : i;
+ all_match &= slab.range_by_band[j].contains(in_p[i], dt_list[i]);
+ }
+ mask_out[pix_idx] |= all_match;
}
- } else {
- // pixel is NDV only if all bands are NDV
- for(size_t i=0; i<nsamps; i++) {
- if(!band_mask[i]) total_mask[i] = 0;
+
+ if(invert) {
+ mask_out[pix_idx] = mask_out[pix_idx] ? 0 : 1;
+ }
+
+ for(size_t band_idx=0; band_idx<bands.size(); band_idx++) {
+ in_p[band_idx] += dt_sizes[band_idx];
}
}
}
-template void NdvDef::arrayCheckNdv<uint8_t>(
- size_t band, const uint8_t *in_data,
- uint8_t *mask_out, size_t nsamps
-) const;
-
-template void NdvDef::arrayCheckNdv<double>(
- size_t band, const double *in_data,
- uint8_t *mask_out, size_t nsamps
-) const;
+void NdvDef::getNdvMask(
+ const void *band, GDALDataType dt,
+ uint8_t *mask_out, size_t num_pixels
+) const {
+ std::vector<const void *> bands;
+ bands.push_back(band);
+ std::vector<GDALDataType> dt_list;
+ dt_list.push_back(dt);
+ getNdvMask(bands, dt_list, mask_out, num_pixels);
+}
} // namespace dangdal
diff --git a/src/ndv.h b/src/ndv.h
index 881e3bd..d5aedbf 100644
--- a/src/ndv.h
+++ b/src/ndv.h
@@ -33,18 +33,33 @@ This code was developed by Dan Stahlke for the Geographic Information Network of
#include <string>
#include <utility>
#include <vector>
+#include <complex>
+
+#include <boost/foreach.hpp>
+
+#include <gdal.h>
+
+#include "datatype_conversion.h"
namespace dangdal {
struct NdvInterval : std::pair<double, double> {
NdvInterval() { }
- NdvInterval(double _min, double _max) :
+ NdvInterval(double _min, double _max) :
std::pair<double, double>(_min, _max) { }
explicit NdvInterval(const std::string &s);
- bool contains(double v) const {
+ template <typename T>
+ bool contains(T v) const {
return v >= first && v <= second;
}
+
+ template <typename T>
+ bool contains(std::complex<T> v) const {
+ return v.real() >= first && v.real() <= second;
+ }
+
+ bool contains(const void *, GDALDataType) const;
};
struct NdvSlab {
@@ -64,18 +79,48 @@ public:
bool empty() const { return slabs.empty(); }
bool isInvert() const { return invert; }
- template<class T>
- void arrayCheckNdv(
- size_t band, const T *in_data,
- uint8_t *mask_out, size_t nsamps
+ void getNdvMask(
+ const void *band, GDALDataType dt,
+ uint8_t *mask_out, size_t num_pixels
) const;
- void aggregateMask(
- uint8_t *total_mask,
- const uint8_t *band_mask,
- size_t nsamps
+ void getNdvMask(
+ const std::vector<const void *> &bands,
+ const std::vector<GDALDataType> &dt_list,
+ uint8_t *mask_out, size_t num_pixels
) const;
+ template <typename T>
+ void getNdvMask(
+ const std::vector<std::vector<T> > &bands,
+ uint8_t *mask_out, size_t num_pixels
+ ) const {
+ GDALDataType dt = GetGDALDataTypeFor<T>::t;
+
+ std::vector<const void *> band_p;
+ std::vector<GDALDataType> dt_list;
+ BOOST_FOREACH(const std::vector<T> &v, bands) {
+ band_p.push_back(reinterpret_cast<const void *>(&v[0]));
+ dt_list.push_back(dt);
+ }
+ getNdvMask(band_p, dt_list, mask_out, num_pixels);
+ }
+
+ // For this one, it doesn't matter what T is, it is the GDALDataType that determines how
+ // the data will be read from memory (i.e. we cast &bands[i][0] to (void *)).
+ template <typename T>
+ void getNdvMask(
+ const std::vector<std::vector<T> > &bands,
+ const std::vector<GDALDataType> &dt_list,
+ uint8_t *mask_out, size_t num_pixels
+ ) const {
+ std::vector<const void *> band_p;
+ BOOST_FOREACH(const std::vector<T> &v, bands) {
+ band_p.push_back(reinterpret_cast<const void *>(&v[0]));
+ }
+ getNdvMask(band_p, dt_list, mask_out, num_pixels);
+ }
+
bool invert;
std::vector<NdvSlab> slabs;
};
diff --git a/src/polygon.cc b/src/polygon.cc
index 250b496..67e94cc 100644
--- a/src/polygon.cc
+++ b/src/polygon.cc
@@ -778,7 +778,10 @@ void Mpoly::xy2ll_with_interp(const GeoRef &georef, double toler) {
}
if(need_midpt) {
- if(num_consec++ > 10) fatal_error("convergence error in mpoly_xy2ll_with_interp");
+ //if(num_consec++ > 10) printf("n=%d\n", num_consec);
+ // If you get this error (probably for lines near the poles), increase this
+ // number and then email me.
+ if(num_consec++ > 100) fatal_error("convergence error in mpoly_xy2ll_with_interp");
if(VERBOSE) {
printf(" xy=[%lf,%lf]:[%lf,%lf]\n", xy1.x, xy1.y, xy2.x, xy2.y);
@@ -847,4 +850,62 @@ Mpoly mpoly_from_wktfile(const std::string &fn) {
return ogr_to_mpoly(geom);
}
+void Ring::debug_dump_binary(FILE *fh) const {
+ uint32_t magic = 0x676e6952; // "Ring"
+ fwrite(&magic, sizeof(magic), 1, fh);
+ uint8_t is_hole_byte = is_hole;
+ fwrite(&is_hole_byte, sizeof(is_hole_byte), 1, fh);
+ fwrite(&parent_id, sizeof(parent_id), 1, fh);
+ size_t npts = pts.size();
+ fwrite(&npts, sizeof(npts), 1, fh);
+ fwrite(&pts[0], sizeof(Vertex), npts, fh);
+}
+
+Ring Ring::debug_load_binary(FILE *fh) {
+ Ring ret;
+
+ uint32_t magic;
+ fread(&magic, sizeof(magic), 1, fh);
+ if(magic != 0x676e6952) {
+ fatal_error("couldn't read ring: bad magic number");
+ }
+ uint8_t is_hole_byte;
+ fread(&is_hole_byte, sizeof(is_hole_byte), 1, fh);
+ ret.is_hole = is_hole_byte;
+ fread(&ret.parent_id, sizeof(ret.parent_id), 1, fh);
+ size_t npts;
+ fread(&npts, sizeof(npts), 1, fh);
+ ret.pts.resize(npts);
+ fread(&ret.pts[0], sizeof(Vertex), npts, fh);
+
+ return ret;
+}
+
+void Mpoly::debug_dump_binary(FILE *fh) const {
+ uint32_t magic = 0x796c704d; // "Mply"
+ fwrite(&magic, sizeof(magic), 1, fh);
+ size_t nrings = rings.size();
+ fwrite(&nrings, sizeof(nrings), 1, fh);
+ for(size_t i=0; i<nrings; i++) {
+ rings[i].debug_dump_binary(fh);
+ }
+}
+
+Mpoly Mpoly::debug_load_binary(FILE *fh) {
+ Mpoly ret;
+
+ uint32_t magic;
+ fread(&magic, sizeof(magic), 1, fh);
+ if(magic != 0x796c704d) {
+ fatal_error("couldn't read mply: bad magic number");
+ }
+ size_t nrings;
+ fread(&nrings, sizeof(nrings), 1, fh);
+ for(size_t i=0; i<nrings; i++) {
+ ret.rings.push_back(Ring::debug_load_binary(fh));
+ }
+
+ return ret;
+}
+
} // namespace dangdal
diff --git a/src/polygon.h b/src/polygon.h
index b29be84..8c983bd 100644
--- a/src/polygon.h
+++ b/src/polygon.h
@@ -32,6 +32,8 @@ This code was developed by Dan Stahlke for the Geographic Information Network of
#include <vector>
#include <string>
#include <algorithm>
+#include <utility>
+#include <cassert>
#include <ogr_api.h>
@@ -63,6 +65,21 @@ public:
empty(true)
{ }
+ Bbox(
+ double _min_x, double _max_x,
+ double _min_y, double _max_y
+ ) :
+ min_x(_min_x), max_x(_max_x),
+ min_y(_min_y), max_y(_max_y),
+ empty(false)
+ { }
+
+ bool contains(const Vertex &v) const {
+ return !empty &&
+ v.x >= min_x && v.x <= max_x &&
+ v.y >= min_y && v.y <= max_y;
+ }
+
void expand(const Vertex &v) {
if(empty) {
empty = false;
@@ -78,6 +95,9 @@ public:
void expand(const Bbox &bb);
+ double width() { return max_x - min_x; }
+ double height() { return max_y - min_y; }
+
double min_x, max_x, min_y, max_y;
bool empty;
};
@@ -93,6 +113,140 @@ static inline bool is_disjoint(const Bbox &bb1, const Bbox &bb2) {
bb2.min_y > bb1.max_y;
}
+// BboxBinarySpacePartition helps you quickly find which of a list of items (that have bounding
+// boxes) intersect a given bounding box. Think of it as a std::map whose keys are bounding
+// boxes. Since several items may intersect a given query box, query returns a list of
+// matches.
+//
+// When many object have overlapping bbox, this tree may not be so efficient and maybe
+// something better should be cooked up if it turns out to be slow.
+template <typename T>
+class BboxBinarySpacePartition {
+public:
+ BboxBinarySpacePartition(
+ std::vector<std::pair<Bbox, T> > items,
+ size_t max_leaf_size = 20,
+ bool axis = 0
+ ) :
+ left(NULL),
+ right(NULL)
+ {
+ leaf_items = items;
+
+ // Compute bbox containing all items.
+ for(size_t i=0; i<items.size(); i++) {
+ node_bbox.expand(items[i].first);
+ }
+
+ if(items.size() > max_leaf_size) {
+ subdivide(max_leaf_size, axis);
+ }
+ }
+
+ ~BboxBinarySpacePartition() {
+ delete(left);
+ delete(right);
+ }
+
+ std::vector<T> get_intersecting_items(Bbox needle) const {
+ std::vector<T> ret;
+ append_intersecting_items(ret, needle);
+ return ret;
+ }
+
+private:
+ void append_intersecting_items(
+ std::vector<T> &out,
+ Bbox needle
+ ) const {
+ if(left) {
+ assert(right);
+ if(!is_disjoint(needle, left->node_bbox)) {
+ left ->append_intersecting_items(out, needle);
+ }
+ if(!is_disjoint(needle, right->node_bbox)) {
+ right->append_intersecting_items(out, needle);
+ }
+ } else {
+ for(size_t i=0; i<leaf_items.size(); i++) {
+ if(!is_disjoint(needle, leaf_items[i].first)) {
+ out.push_back(leaf_items[i].second);
+ }
+ }
+ }
+ }
+
+ void subdivide(size_t max_leaf_size, bool axis) {
+ std::vector<double> midpts(leaf_items.size());
+ double cmin = 0, cmax = 0, cavg = 0;
+ size_t num_nonempty = 0;
+ for(size_t i=0; i<leaf_items.size(); i++) {
+ const Bbox &ibb = leaf_items[i].first;
+ if(ibb.empty) continue;
+ num_nonempty++;
+ node_bbox.expand(ibb);
+ double item_mid = axis ?
+ (ibb.min_y + ibb.max_y) / 2.0 :
+ (ibb.min_x + ibb.max_x) / 2.0;
+ if(i == 0) cmin = cmax = item_mid;
+ cmin = std::min(cmin, item_mid);
+ cmax = std::max(cmax, item_mid);
+ cavg += item_mid;
+ }
+ cavg /= num_nonempty;
+
+ // Subdivide space into two halves. Note: this is used to decide which branch of the
+ // tree each item goes into. However, each branch will then compute the exact Bbox
+ // that contains each of its items (it won't use these ones we compute here).
+ Bbox left_bbox, right_bbox;
+ if(axis) {
+ left_bbox = Bbox(
+ node_bbox.min_x, node_bbox.max_x,
+ cmin, cavg);
+ right_bbox = Bbox(
+ node_bbox.min_x, node_bbox.max_x,
+ cavg, cmax);
+ } else {
+ left_bbox = Bbox(
+ cmin, cavg,
+ node_bbox.min_y, node_bbox.max_y);
+ right_bbox = Bbox(
+ cavg, cmax,
+ node_bbox.min_y, node_bbox.max_y);
+ }
+
+ // Find which items go in each box.
+ //
+ // Note: it would be possible to just move items to the left and right side of the
+ // original array like with quicksort, and never have to allocate all of these little
+ // vectors for every node.
+ std::vector<std::pair<Bbox, T> > left_items, right_items;
+ for(size_t i=0; i<leaf_items.size(); i++) {
+ const Bbox &ibb = leaf_items[i].first;
+ Vertex item_mid = Vertex(
+ (ibb.min_x + ibb.max_x) / 2.0,
+ (ibb.min_y + ibb.max_y) / 2.0);
+ if(left_bbox.contains(item_mid)) {
+ left_items.push_back(leaf_items[i]);
+ } else {
+ right_items.push_back(leaf_items[i]);
+ }
+ }
+
+ // free memory
+ leaf_items.clear();
+
+ left = new BboxBinarySpacePartition<T>(left_items, max_leaf_size, !axis);
+ right = new BboxBinarySpacePartition<T>(right_items, max_leaf_size, !axis);
+ }
+
+private:
+ Bbox node_bbox;
+ BboxBinarySpacePartition<T> *left, *right;
+ // only used by leafs
+ std::vector<std::pair<Bbox, T> > leaf_items;
+};
+
class Ring {
public:
Ring() : is_hole(false), parent_id(-1) { }
@@ -112,6 +266,9 @@ public:
return ret;
}
+ void debug_dump_binary(FILE *fh) const;
+ static Ring debug_load_binary(FILE *fh);
+
std::vector<Vertex> pts;
bool is_hole;
int parent_id;
@@ -131,6 +288,9 @@ public:
void en2xy(const GeoRef &georef);
void xy2ll_with_interp(const GeoRef &georef, double toler);
+ void debug_dump_binary(FILE *fh) const;
+ static Mpoly debug_load_binary(FILE *fh);
+
std::vector<Ring> rings;
};
diff --git a/src/raster_features.cc b/src/raster_features.cc
new file mode 100644
index 0000000..8502c8c
--- /dev/null
+++ b/src/raster_features.cc
@@ -0,0 +1,440 @@
+/*
+Copyright (c) 2013, Regents of the University of Alaska
+
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
+
+ * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
+ * Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
+ * Neither the name of the Geographic Information Network of Alaska nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
+CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+This code was developed by Dan Stahlke for the Geographic Information Network of Alaska.
+*/
+
+
+
+#include <boost/foreach.hpp>
+
+#include "raster_features.h"
+
+namespace dangdal {
+
+FeatureInterpreter::BandInfo::BandInfo() :
+ raw_val_offset(0),
+ raw_val_size(0),
+ is_int(0),
+ is_double(0),
+ is_complex(0),
+ color_table(NULL),
+ ogr_fld(NULL)
+{ }
+
+// GDALCopyWords doesn't allow const input
+int32_t FeatureInterpreter::BandInfo::val_to_int(FeatureRawVal rawval) const {
+ return gdal_scalar_to_int32(&rawval[raw_val_offset], dt);
+}
+
+// GDALCopyWords doesn't allow const input
+double FeatureInterpreter::BandInfo::val_to_double(FeatureRawVal rawval) const {
+ return gdal_scalar_to_double(&rawval[raw_val_offset], dt);
+}
+
+std::string FeatureInterpreter::BandInfo::val_to_str(const FeatureRawVal &rawval) const {
+ if(color_table) {
+ const GDALColorEntry *color = GDALGetColorEntry(color_table,
+ val_to_int(rawval));
+ return str(boost::format("%d, %d, %d, %d") %
+ color->c1 % color->c2 % color->c3 % color->c4);
+ } else if(is_int) {
+ return str(boost::format("%d") % val_to_int(rawval));
+ } else if(is_double) {
+ return str(boost::format("%g") % val_to_double(rawval));
+ } else {
+ std::string ret = "0x";
+ for(size_t i=0; i<raw_val_size; i++) {
+ uint8_t v = rawval[raw_val_offset + i];
+ ret += str(boost::format("%02x") % int(v));
+ }
+ return ret;
+ }
+}
+
+FeatureInterpreter::FeatureInterpreter(GDALDatasetH ds, std::vector<size_t> band_ids) {
+ size_t offset = 0;
+ BOOST_FOREACH(const size_t band_id, band_ids) {
+ BandInfo bi;
+
+ GDALRasterBandH band = GDALGetRasterBand(ds, band_id);
+ if(!band) fatal_error("Could not open band %zd.", band_id);
+
+ if(band_ids.size() == 1) {
+ bi.label = "value";
+ } else {
+ bi.label = str(boost::format("band%d") % band_id);
+ }
+
+ bi.dt = GDALGetRasterDataType(band);
+ bi.raw_val_size = GDALGetDataTypeSize(bi.dt) / 8;
+ bi.raw_val_offset = offset;
+ offset += bi.raw_val_size;
+
+ bi.is_int = bi.is_double = bi.is_complex = false;
+ switch(bi.dt) {
+ case GDT_Byte:
+ case GDT_UInt16:
+ case GDT_Int16:
+ case GDT_UInt32:
+ case GDT_Int32:
+ bi.is_int = true; break;
+ case GDT_Float32:
+ case GDT_Float64:
+ bi.is_double = true; break;
+ case GDT_CInt16:
+ case GDT_CInt32:
+ case GDT_CFloat32:
+ case GDT_CFloat64:
+ bi.is_complex = true; break;
+ case GDT_Unknown:
+ default:
+ break; // no-op
+ }
+
+ switch(bi.dt) {
+ case GDT_Byte:
+ case GDT_UInt16:
+ case GDT_Int16:
+ case GDT_UInt32:
+ case GDT_Int32:
+ bi.ogr_ft = OFTInteger; break;
+ case GDT_Float32:
+ case GDT_Float64:
+ bi.ogr_ft = OFTReal; break;
+ case GDT_CInt16:
+ case GDT_CInt32:
+ case GDT_CFloat32:
+ case GDT_CFloat64:
+ // FIXME - currently complex goes to string, but a pair of numeric fields
+ // would be better.
+ bi.ogr_ft = OFTString; break;
+ case GDT_Unknown:
+ default:
+ bi.ogr_ft = OFTString; break;
+ }
+ bi.ogr_fld = OGR_Fld_Create(bi.label.c_str(), bi.ogr_ft);
+
+ if(GDALGetRasterColorInterpretation(band) == GCI_PaletteIndex) {
+ bi.color_table = GDALGetRasterColorTable(band);
+ for(int i=0; i<gdal_color_table_size; i++) {
+ std::string label;
+ if(band_ids.size() > 1) {
+ label = str(boost::format("band%d.") % band_id);
+ }
+ label += str(boost::format("c%d") % (i+1));
+
+ bi.ogr_pal_fld[i].first = label;
+ bi.ogr_pal_fld[i].second = OGR_Fld_Create(label.c_str(), OFTInteger);
+ }
+ } else {
+ bi.color_table = NULL;
+ }
+
+ band_info_list.push_back(bi);
+ }
+}
+
+std::string FeatureInterpreter::pixel_to_string(const FeatureRawVal &rawval) const {
+ std::string ret;
+ for(size_t i=0; i<band_info_list.size(); i++) {
+ if(i) ret += ", ";
+ ret += band_info_list[i].val_to_str(rawval);
+ }
+ return ret;
+}
+
+void FeatureInterpreter::create_ogr_fields(OGRLayerH ogr_layer) const {
+ BOOST_FOREACH(const BandInfo &bi, band_info_list) {
+ OGR_L_CreateField(ogr_layer, bi.ogr_fld, TRUE);
+ if(bi.color_table) {
+ typedef std::map<std::string, OGRFieldDefnH>::value_type fld_pair_t;
+ BOOST_FOREACH(const fld_pair_t &f, bi.ogr_pal_fld) {
+ OGR_L_CreateField(ogr_layer, f.second, TRUE);
+ }
+ }
+ }
+}
+
+void FeatureInterpreter::set_ogr_fields(OGRLayerH ogr_layer, OGRFeatureH ogr_feat, const FeatureRawVal &rawval) const {
+ BOOST_FOREACH(const BandInfo &bi, band_info_list) {
+ int fld_idx = OGR_FD_GetFieldIndex(OGR_L_GetLayerDefn(ogr_layer), bi.label.c_str());
+ // Did you remember to call create_ogr_fields?
+ if(fld_idx < 0) fatal_error("cannot get OGR field %s", bi.label.c_str());
+ switch(bi.ogr_ft) {
+ case OFTInteger:
+ OGR_F_SetFieldInteger(ogr_feat, fld_idx, bi.val_to_int(rawval));
+ break;
+ case OFTReal:
+ OGR_F_SetFieldDouble(ogr_feat, fld_idx, bi.val_to_double(rawval));
+ break;
+ case OFTString:
+ OGR_F_SetFieldString(ogr_feat, fld_idx, bi.val_to_str(rawval).c_str());
+ break;
+ default:
+ // This should not happen, since ogr_ft is only set by this class.
+ fatal_error("not implemented");
+ }
+
+ if(bi.color_table) {
+ const GDALColorEntry *color = GDALGetColorEntry(bi.color_table, bi.val_to_int(rawval));
+ for(int i=0; i<gdal_color_table_size; i++) {
+ const std::pair<std::string, OGRFieldDefnH> &f = bi.ogr_pal_fld[i];
+ int fld_idx = OGR_FD_GetFieldIndex(OGR_L_GetLayerDefn(ogr_layer), f.first.c_str());
+ if(fld_idx < 0) fatal_error("cannot get OGR field %s", f.first.c_str());
+ int val =
+ i==0 ? color->c1 :
+ i==1 ? color->c2 :
+ i==2 ? color->c3 : color->c4;
+ OGR_F_SetFieldInteger(ogr_feat, fld_idx, val);
+ }
+ }
+ }
+}
+
+FeatureBitmap::FeatureBitmap(const size_t _w, const size_t _h, const size_t _raw_vals_size) :
+ w(_w), h(_h),
+ raw_vals_size(_raw_vals_size),
+ raster(w, h)
+{ }
+
+FeatureBitmap::Index FeatureBitmap::get_index(const FeatureRawVal &pixel) {
+ std::map<FeatureRawVal, Index>::iterator it = table.find(pixel);
+ if(it == table.end()) {
+ if(table.size() >= std::numeric_limits<FeatureBitmap::Index>::max()) {
+ fatal_error("Input had too many feature values (max is %zd)",
+ size_t(std::numeric_limits<FeatureBitmap::Index>::max()));
+ }
+ FeatureBitmap::Index v = table.size();
+ table[pixel] = v;
+ return v;
+ } else {
+ return it->second;
+ }
+}
+
+void FeatureBitmap::dump_feature_table() const {
+ typedef std::map<FeatureRawVal, Index>::value_type table_pair_t;
+ BOOST_FOREACH(const table_pair_t &f, table) {
+ printf("feature %d:", f.second);
+ for(size_t i=0; i<f.first.size(); i++) {
+ if(i) printf(",");
+ printf(" %d", f.first[i]);
+ }
+ printf("\n");
+ }
+}
+
+BitGrid FeatureBitmap::get_mask_for_feature(FeatureBitmap::Index wanted) const {
+ BitGrid mask(w, h);
+
+ for(size_t y=0; y<h; y++) {
+ for(size_t x=0; x<w; x++) {
+ mask.set(x, y, raster(x, y) == wanted);
+ }
+ }
+
+ return mask;
+}
+
+FeatureBitmap *FeatureBitmap::from_raster(
+ GDALDatasetH ds, std::vector<size_t> band_ids, const NdvDef &ndv_def, DebugPlot *dbuf
+) {
+ assert(!band_ids.empty());
+
+ size_t w = GDALGetRasterXSize(ds);
+ size_t h = GDALGetRasterYSize(ds);
+ size_t band_count = GDALGetRasterCount(ds);
+ if(VERBOSE) printf("input is %zd x %zd x %zd\n", w, h, band_count);
+
+ std::vector<GDALRasterBandH> bands;
+ BOOST_FOREACH(const size_t band_id, band_ids) {
+ if(VERBOSE) printf("opening band %zd\n", band_id);
+ GDALRasterBandH band = GDALGetRasterBand(ds, band_id);
+ if(!band) fatal_error("Could not open band %zd.", band_id);
+ bands.push_back(band);
+ }
+
+ int blocksize_x_int, blocksize_y_int;
+ GDALGetBlockSize(bands[0], &blocksize_x_int, &blocksize_y_int);
+ // Out of laziness, I am hoping that images always have the same block size for each band.
+ BOOST_FOREACH(const GDALRasterBandH band, bands) {
+ int bx, by;
+ GDALGetBlockSize(band, &bx, &by);
+ if(bx != blocksize_x_int || by != blocksize_y_int) {
+ fatal_error(
+ "Bands have different block sizes. Not currently implemented. Please contact developer");
+ }
+ }
+ size_t blocksize_x = blocksize_x_int;
+ size_t blocksize_y = blocksize_y_int;
+ size_t blocksize_xy = blocksize_x * blocksize_y;
+
+ std::vector<GDALDataType> datatypes;
+ std::vector<size_t> dt_sizes;
+ size_t dt_total_size = 0;
+ if(VERBOSE >= 2) printf("datatype sizes:");
+ BOOST_FOREACH(const GDALRasterBandH band, bands) {
+ GDALDataType dt = GDALGetRasterDataType(band);
+ datatypes.push_back(dt);
+ size_t s = GDALGetDataTypeSize(dt) / 8;
+ dt_sizes.push_back(s);
+ dt_total_size += s;
+ if(VERBOSE >= 2) printf(" %zd", s);
+ }
+ if(VERBOSE >= 2) printf("\n");
+
+ std::vector<std::vector<uint8_t> > band_buf(bands.size());
+ for(size_t i=0; i<bands.size(); i++) {
+ size_t num_bytes = blocksize_xy * dt_sizes[i];
+ band_buf[i].resize(num_bytes);
+ }
+
+ std::vector<uint8_t> ndv_mask(blocksize_xy);
+
+ FeatureBitmap *fbm = new FeatureBitmap(w, h, dt_total_size);
+
+ size_t num_valid = 0;
+ size_t num_ndv = 0;
+
+ printf("Reading input...\n");
+
+ size_t num_blocks_x = (w + blocksize_x - 1) / blocksize_x;
+ size_t num_blocks_y = (h + blocksize_y - 1) / blocksize_y;
+ for(size_t block_y=0; block_y<num_blocks_y; block_y++) {
+ size_t boff_y = blocksize_y * block_y;
+ size_t bsize_y = blocksize_y;
+ if(bsize_y + boff_y > h) bsize_y = h - boff_y;
+ for(size_t block_x=0; block_x<num_blocks_x; block_x++) {
+ size_t boff_x = blocksize_x * block_x;
+ size_t bsize_x = blocksize_x;
+ if(bsize_x + boff_x > w) bsize_x = w - boff_x;
+
+ double progress =
+ double(
+ boff_y * w +
+ boff_x * bsize_y
+ ) / (w * h);
+ GDALTermProgress(progress, NULL, NULL);
+
+ for(size_t band_idx=0; band_idx<bands.size(); band_idx++) {
+ GDALReadBlock(bands[band_idx], block_x, block_y, &band_buf[band_idx][0]);
+ }
+ if(!ndv_def.empty()) {
+ ndv_def.getNdvMask(band_buf, datatypes, &ndv_mask[0], blocksize_xy);
+ }
+
+ FeatureRawVal pixel;
+ pixel.resize(dt_total_size);
+
+ for(size_t sub_y=0; sub_y<bsize_y; sub_y++) {
+ size_t y = sub_y + boff_y;
+ bool is_dbuf_stride_y = dbuf && ((y % dbuf->stride_y) == 0);
+
+ for(size_t sub_x=0; sub_x<bsize_x; sub_x++) {
+ size_t x = sub_x + boff_x;
+ bool is_dbuf_stride = is_dbuf_stride_y && ((sub_x % dbuf->stride_x) == 0);
+
+ size_t in_idx = blocksize_x*sub_y + sub_x;
+
+ if(ndv_mask[in_idx]) {
+ num_ndv++;
+
+ if(is_dbuf_stride) {
+ dbuf->plotPoint(x, y, 0, 0, 0);
+ }
+ } else {
+ num_valid++;
+
+ size_t j = 0;
+ for(size_t band_id=0; band_id<bands.size(); band_id++) {
+ uint8_t *p = &band_buf[band_id][in_idx*dt_sizes[band_id]];
+ for(size_t i=0; i<dt_sizes[band_id]; i++) {
+ pixel[j++] = *(p++);
+ }
+ }
+ assert(j == dt_total_size);
+
+ FeatureBitmap::Index index_val = fbm->get_index(pixel);
+ fbm->raster(x, y) = index_val;
+
+ if(is_dbuf_stride) {
+ size_t x = sub_x + boff_x;
+ // assign a random palette for the debug report
+ uint8_t r = index_val * 100 + 100;
+ uint8_t g = index_val * 173 + 36;
+ uint8_t b = index_val * 47 + 202;
+ dbuf->plotPoint(x, y, r, g, b);
+ }
+ }
+ }
+ }
+ }
+ }
+
+ GDALTermProgress(1, NULL, NULL);
+
+ printf("Found %zd valid and %zd NDV pixels.\n", num_valid, num_ndv);
+
+ return fbm;
+}
+
+} // namespace dangdal
+
+//using namespace dangdal;
+//void usage(const std::string &cmdname) {}
+//int main(int argc, char **argv) {
+// if(argc != 2) fatal_error("give filename");
+//
+// VERBOSE = 2;
+//
+// GDALAllRegister();
+//
+// std::string input_raster_fn(argv[1]);
+// GDALDatasetH ds = GDALOpen(input_raster_fn.c_str(), GA_ReadOnly);
+// if(!ds) fatal_error("open failed");
+// size_t w = GDALGetRasterXSize(ds);
+// size_t h = GDALGetRasterYSize(ds);
+//
+// std::vector<size_t> band_ids;
+// size_t nbands = GDALGetRasterCount(ds);
+// for(size_t i=0; i<nbands; i++) band_ids.push_back(i+1);
+//
+// DebugPlot *dbuf = new DebugPlot(w, h, PLOT_CONTOURS);
+//
+// NdvDef ndv_def = NdvDef(ds, band_ids);
+//
+// FeatureBitmap *fbm = FeatureBitmap::from_raster(ds, band_ids, ndv_def, dbuf);
+// FeatureInterpreter interp(ds, band_ids);
+//
+// dbuf->writePlot("zz.ppm");
+//
+// typedef std::map<FeatureRawVal, FeatureBitmap::Index>::value_type feature_pair_t;
+// BOOST_FOREACH(const feature_pair_t &f, fbm->feature_table()) {
+// printf("feature %d: %s\n", f.second, interp.pixel_to_string(f.first).c_str());
+// }
+//
+// GDALClose(ds);
+//
+// return 0;
+//}
diff --git a/src/raster_features.h b/src/raster_features.h
new file mode 100644
index 0000000..c42bd96
--- /dev/null
+++ b/src/raster_features.h
@@ -0,0 +1,118 @@
+/*
+Copyright (c) 2013, Regents of the University of Alaska
+
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
+
+ * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
+ * Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
+ * Neither the name of the Geographic Information Network of Alaska nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
+CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+This code was developed by Dan Stahlke for the Geographic Information Network of Alaska.
+*/
+
+
+
+// This module supports the gdal_trace_outline -classify option. A raster is read and pixel
+// values are converted to unique indices (kind of like a palette, except for arbitrary data
+// types). With this indexed array, it is possible to pull out a bitmask for each individual
+// pixel value. Also, the pixel values can be formatted as strings or as fields in an OGR
+// file.
+
+#include <vector>
+#include <map>
+#include <utility>
+#include <limits>
+#include <string>
+
+#include <boost/format.hpp>
+
+#include <ogr_api.h>
+
+#include "common.h"
+#include "mask.h"
+#include "debugplot.h"
+#include "datatype_conversion.h"
+
+namespace dangdal {
+
+// Represents a raw pixel value. For convenience, it inherits directly from std::vector so as
+// to inherit the comparison operators.
+struct FeatureRawVal : public std::vector<uint8_t> { };
+
+// Converts pixel values to strings or OGR fields.
+struct FeatureInterpreter {
+private:
+ static const int gdal_color_table_size = 4;
+
+ struct BandInfo {
+ BandInfo();
+
+ // GDALCopyWords doesn't allow const input
+ int32_t val_to_int(FeatureRawVal rawval) const;
+ // GDALCopyWords doesn't allow const input
+ double val_to_double(FeatureRawVal rawval) const;
+
+ std::string val_to_str(const FeatureRawVal &rawval) const;
+
+ std::string label;
+ size_t raw_val_offset;
+ size_t raw_val_size;
+ GDALDataType dt;
+ bool is_int, is_double, is_complex;
+ GDALColorTableH color_table;
+ OGRFieldType ogr_ft;
+ OGRFieldDefnH ogr_fld;
+ // For paletted bands, this gives the label and field for each color component.
+ std::pair<std::string, OGRFieldDefnH> ogr_pal_fld[gdal_color_table_size];
+ };
+
+public:
+ FeatureInterpreter(GDALDatasetH ds, std::vector<size_t> band_ids);
+ std::string pixel_to_string(const FeatureRawVal &rawval) const;
+ void create_ogr_fields(OGRLayerH ogr_layer) const;
+ void set_ogr_fields(OGRLayerH ogr_layer, OGRFeatureH ogr_feat, const FeatureRawVal &rawval) const;
+
+private:
+ std::vector<BandInfo> band_info_list;
+};
+
+// A bitmap of features. Features are stored as type Index in a bitmap, and can be mapped to
+// FeatureRawVal.
+struct FeatureBitmap {
+ typedef uint16_t Index;
+
+ FeatureBitmap(const size_t _w, const size_t _h, const size_t _raw_vals_size);
+
+ static FeatureBitmap *from_raster(
+ GDALDatasetH ds, std::vector<size_t> band_ids, const NdvDef &ndv_def, DebugPlot *dbuf);
+
+ const std::map<FeatureRawVal, Index> &feature_table() const {
+ return table;
+ }
+
+ Index get_index(const FeatureRawVal &pixel);
+ void dump_feature_table() const;
+ BitGrid get_mask_for_feature(Index wanted) const;
+
+private:
+ const size_t w, h;
+ const size_t raw_vals_size;
+ GridArray<Index> raster;
+ std::map<FeatureRawVal, Index> table;
+};
+
+} // namespace dangdal
diff --git a/src/tests/.gitignore b/src/tests/.gitignore
index edd53c6..d2259e0 100644
--- a/src/tests/.gitignore
+++ b/src/tests/.gitignore
@@ -1,2 +1,6 @@
# testcase output
out_*
+
+# generated as needed
+gradient[123].tif
+has_nan.tif
diff --git a/src/tests/good_test1_1.ppm b/src/tests/good_test1_1.ppm
new file mode 100644
index 0000000..d99bda4
Binary files /dev/null and b/src/tests/good_test1_1.ppm differ
diff --git a/src/tests/good_test2_2.ppm b/src/tests/good_test1_2.ppm
similarity index 96%
copy from src/tests/good_test2_2.ppm
copy to src/tests/good_test1_2.ppm
index e308f38..d0263e9 100644
Binary files a/src/tests/good_test2_2.ppm and b/src/tests/good_test1_2.ppm differ
diff --git a/src/tests/good_test1_3.ppm b/src/tests/good_test1_3.ppm
new file mode 100644
index 0000000..91f213e
Binary files /dev/null and b/src/tests/good_test1_3.ppm differ
diff --git a/src/tests/good_test1_3_classify.dbf b/src/tests/good_test1_3_classify.dbf
new file mode 100644
index 0000000..bbbe528
Binary files /dev/null and b/src/tests/good_test1_3_classify.dbf differ
diff --git a/src/tests/good_test1_3_classify.shp b/src/tests/good_test1_3_classify.shp
new file mode 100644
index 0000000..43e2f68
Binary files /dev/null and b/src/tests/good_test1_3_classify.shp differ
diff --git a/src/tests/good_test1_3_classify.shx b/src/tests/good_test1_3_classify.shx
new file mode 100644
index 0000000..3450a25
Binary files /dev/null and b/src/tests/good_test1_3_classify.shx differ
diff --git a/src/tests/good_test1_3_classify_pal.dbf b/src/tests/good_test1_3_classify_pal.dbf
new file mode 100644
index 0000000..a5648a6
Binary files /dev/null and b/src/tests/good_test1_3_classify_pal.dbf differ
diff --git a/src/tests/good_test1_3_classify_pal.shp b/src/tests/good_test1_3_classify_pal.shp
new file mode 100644
index 0000000..487734c
Binary files /dev/null and b/src/tests/good_test1_3_classify_pal.shp differ
diff --git a/src/tests/good_test1_3_classify_pal.shx b/src/tests/good_test1_3_classify_pal.shx
new file mode 100644
index 0000000..065b75b
Binary files /dev/null and b/src/tests/good_test1_3_classify_pal.shx differ
diff --git a/src/tests/good_test1_4-rect.ppm b/src/tests/good_test1_4-rect.ppm
index a9fb147..d6c8992 100644
Binary files a/src/tests/good_test1_4-rect.ppm and b/src/tests/good_test1_4-rect.ppm differ
diff --git a/src/tests/good_test2_4.ppm b/src/tests/good_test1_4.ppm
similarity index 71%
copy from src/tests/good_test2_4.ppm
copy to src/tests/good_test1_4.ppm
index a5e04c1..6ad27b0 100644
Binary files a/src/tests/good_test2_4.ppm and b/src/tests/good_test1_4.ppm differ
diff --git a/src/tests/good_test2_5.ppm b/src/tests/good_test1_5.ppm
similarity index 65%
copy from src/tests/good_test2_5.ppm
copy to src/tests/good_test1_5.ppm
index 09b8d08..78a6bef 100644
Binary files a/src/tests/good_test2_5.ppm and b/src/tests/good_test1_5.ppm differ
diff --git a/src/tests/good_test1_double.dbf b/src/tests/good_test1_double.dbf
new file mode 100644
index 0000000..8c95ac5
Binary files /dev/null and b/src/tests/good_test1_double.dbf differ
diff --git a/src/tests/good_test1_double.shp b/src/tests/good_test1_double.shp
new file mode 100644
index 0000000..b2c55e2
Binary files /dev/null and b/src/tests/good_test1_double.shp differ
diff --git a/src/tests/good_test1_double.shx b/src/tests/good_test1_double.shx
new file mode 100644
index 0000000..c244a26
Binary files /dev/null and b/src/tests/good_test1_double.shx differ
diff --git a/src/tests/good_test1_double.wkt b/src/tests/good_test1_double.wkt
new file mode 100644
index 0000000..9234841
--- /dev/null
+++ b/src/tests/good_test1_double.wkt
@@ -0,0 +1,13 @@
+MULTIPOLYGON (((157 87,166 87,166 88,169 88,169 89,171 89,171 90,180 90,180 91,185 91,185 92,188 92,188 93,199 93,199 101,198 101,198 103,199 103,199 104,198 104,198 106,203 106,203 107,204 107,204 109,205 109,205 113,208 113,208 112,211 112,211 111,218 111,218 113,219 113,219 114,220 114,220 115,222 115,222 116,224 116,224 117,230 117,230 118,241 118,241 119,249 119,249 120,260 120,260 119,261 119,261 120,263 120,263 119,283 119,283 118,285 118,285 119,284 119,284 120,283 120,283 121,26 [...]
+MULTIPOLYGON (((160 97,164 97,164 98,165 98,165 99,166 99,166 100,167 100,167 101,171 101,171 102,178 102,178 103,182 103,182 105,181 105,181 106,182 106,182 107,183 107,183 108,184 108,184 109,188 109,188 110,190 110,190 111,192 111,192 112,193 112,193 113,195 113,195 114,196 114,196 115,197 115,197 116,199 116,199 117,200 117,200 118,204 118,204 119,206 119,206 120,209 120,209 119,211 119,211 118,215 118,215 119,216 119,216 121,217 121,217 123,218 123,218 124,220 124,220 125,222 125,22 [...]
+MULTIPOLYGON (((189 120,191 120,191 121,193 121,193 122,195 122,195 123,197 123,197 124,198 124,198 126,200 126,200 127,202 127,202 128,215 128,215 129,216 129,216 130,214 130,214 131,212 131,212 132,211 132,211 133,209 133,209 135,208 135,208 140,207 140,207 141,206 141,206 142,205 142,205 143,204 143,204 144,203 144,203 145,197 145,197 146,188 146,188 147,185 147,185 148,182 148,182 149,180 149,180 150,176 150,176 151,174 151,174 152,170 152,170 153,166 153,166 154,165 154,165 155,164 [...]
+MULTIPOLYGON (((149 159,150 159,150 160,149 160,149 159)),((414 227,421 227,421 228,424 228,424 229,423 229,423 230,419 230,419 231,417 231,417 232,416 232,416 233,415 233,415 234,414 234,414 233,413 233,413 229,414 229,414 227)),((410 246,412 246,412 248,413 248,413 250,412 250,412 251,409 251,409 250,408 250,408 248,410 248,410 246)),((390 262,391 262,391 263,394 263,394 264,395 264,395 265,394 265,394 266,393 266,393 267,392 267,392 268,389 268,389 269,388 269,388 272,378 272,378 273, [...]
+POLYGON ((331 269,340 269,340 270,351 270,351 271,355 271,355 272,357 272,357 273,359 273,359 275,356 275,356 276,332 276,332 277,329 277,329 282,330 282,330 283,331 283,331 286,330 286,330 291,331 291,331 293,330 293,330 294,328 294,328 295,327 295,327 296,325 296,325 297,322 297,322 298,319 298,319 297,316 297,316 296,315 296,315 294,314 294,314 293,312 293,312 292,306 292,306 291,304 291,304 289,305 289,305 288,306 288,306 286,303 286,303 285,297 285,297 284,289 284,289 285,284 285,28 [...]
+MULTIPOLYGON (((328 271,340 271,340 272,342 272,342 273,340 273,340 274,330 274,330 273,326 273,326 272,328 272,328 271)),((313 274,320 274,320 275,321 275,321 277,322 277,322 279,323 279,323 281,324 281,324 284,323 284,323 285,320 285,320 284,318 284,318 283,316 283,316 282,314 282,314 281,311 281,311 280,302 280,302 279,301 279,301 278,303 278,303 277,306 277,306 276,308 276,308 275,313 275,313 274)))
+MULTIPOLYGON (((0 0,8 0,8 1,4 1,4 2,3 2,3 1,2 1,2 2,0 2,0 0)),((486 349,488 349,488 350,489 350,489 351,492 351,492 350,493 350,493 351,494 351,494 352,496 352,496 351,499 351,499 350,501 350,501 351,503 351,503 409,475 409,475 408,481 408,481 407,486 407,486 406,487 406,487 405,489 405,489 402,487 402,487 401,484 401,484 400,483 400,483 399,484 399,484 398,485 398,485 395,486 395,486 391,485 391,485 390,484 390,484 388,482 388,482 383,483 383,483 382,484 382,484 380,483 380,483 379,482 [...]
+MULTIPOLYGON (((8 0,135 0,135 1,133 1,133 2,132 2,132 3,131 3,131 4,130 4,130 6,129 6,129 7,125 7,125 8,124 8,124 9,123 9,123 10,122 10,122 11,121 11,121 12,117 12,117 13,115 13,115 14,110 14,110 15,109 15,109 16,107 16,107 17,105 17,105 19,104 19,104 20,105 20,105 23,103 23,103 22,102 22,102 21,98 21,98 20,97 20,97 19,84 19,84 20,83 20,83 21,82 21,82 23,80 23,80 24,78 24,78 25,77 25,77 26,74 26,74 27,71 27,71 28,68 28,68 29,63 29,63 30,62 30,62 32,60 32,60 33,57 33,57 34,55 34,55 35,54 [...]
+MULTIPOLYGON (((135 0,323 0,323 1,321 1,321 2,319 2,319 3,318 3,318 4,317 4,317 5,316 5,316 6,315 6,315 9,316 9,316 11,315 11,315 13,316 13,316 16,314 16,314 15,311 15,311 14,310 14,310 13,308 13,308 12,305 12,305 11,303 11,303 10,297 10,297 9,293 9,293 10,292 10,292 11,289 11,289 12,284 12,284 13,283 13,283 14,281 14,281 15,280 15,280 16,278 16,278 14,277 14,277 13,274 13,274 14,271 14,271 17,267 17,267 18,266 18,266 20,263 20,263 18,262 18,262 17,259 17,259 16,253 16,253 17,252 17,252 [...]
+MULTIPOLYGON (((323 0,367 0,367 2,365 2,365 3,363 3,363 4,362 4,362 5,363 5,363 7,365 7,365 6,368 6,368 5,370 5,370 4,373 4,373 3,378 3,378 4,380 4,380 5,381 5,381 6,382 6,382 7,383 7,383 8,385 8,385 9,387 9,387 10,390 10,390 9,395 9,395 8,393 8,393 7,392 7,392 3,393 3,393 2,392 2,392 1,393 1,393 0,503 0,503 12,502 12,502 13,501 13,501 14,499 14,499 15,498 15,498 16,497 16,497 17,496 17,496 18,493 18,493 19,492 19,492 20,498 20,498 19,502 19,502 18,503 18,503 33,501 33,501 32,491 32,491 [...]
+MULTIPOLYGON (((426 13,429 13,429 14,430 14,430 13,431 13,431 14,438 14,438 15,441 15,441 14,443 14,443 15,444 15,444 16,446 16,446 17,449 17,449 18,455 18,455 19,457 19,457 20,459 20,459 21,460 21,460 20,463 20,463 21,465 21,465 23,466 23,466 27,467 27,467 29,466 29,466 30,465 30,465 32,466 32,466 33,472 33,472 34,473 34,473 35,479 35,479 36,481 36,481 35,485 35,485 34,489 34,489 33,491 33,491 32,501 32,501 33,503 33,503 55,501 55,501 56,500 56,500 57,493 57,493 58,491 58,491 61,492 61, [...]
+MULTIPOLYGON (((369 49,370 49,370 51,371 51,371 52,372 52,372 53,374 53,374 52,375 52,375 51,377 51,377 52,382 52,382 50,385 50,385 51,387 51,387 52,390 52,390 51,391 51,391 50,392 50,392 51,393 51,393 54,395 54,395 52,397 52,397 54,399 54,399 51,398 51,398 49,404 49,404 50,407 50,407 51,412 51,412 52,414 52,414 53,413 53,413 54,412 54,412 55,415 55,415 56,416 56,416 55,417 55,417 54,418 54,418 55,419 55,419 56,420 56,420 57,421 57,421 58,422 58,422 57,423 57,423 56,425 56,425 57,426 57, [...]
+MULTIPOLYGON (((349 76,360 76,360 77,361 77,361 79,368 79,368 78,372 78,372 79,374 79,374 80,381 80,381 79,411 79,411 80,421 80,421 81,429 81,429 82,433 82,433 83,437 83,437 84,438 84,438 85,441 85,441 86,442 86,442 87,441 87,441 88,434 88,434 89,426 89,426 90,425 90,425 91,422 91,422 92,420 92,420 93,415 93,415 94,411 94,411 95,409 95,409 96,406 96,406 97,401 97,401 98,397 98,397 99,394 99,394 100,390 100,390 101,388 101,388 102,386 102,386 103,384 103,384 104,382 104,382 105,381 105,38 [...]
diff --git a/src/tests/good_test1_double_clip.wkt b/src/tests/good_test1_double_clip.wkt
new file mode 100644
index 0000000..a9b88f6
--- /dev/null
+++ b/src/tests/good_test1_double_clip.wkt
@@ -0,0 +1,3 @@
+MULTIPOLYGON (((0 0,323 0,323 1,321 1,321 2,319 2,319 3,318 3,318 4,317 4,317 5,316 5,316 6,315 6,315 9,316 9,316 11,315 11,315 13,316 13,316 16,314 16,314 15,311 15,311 14,310 14,310 13,308 13,308 12,305 12,305 11,303 11,303 10,297 10,297 9,293 9,293 10,292 10,292 11,289 11,289 12,284 12,284 13,283 13,283 14,281 14,281 15,280 15,280 16,278 16,278 14,277 14,277 13,274 13,274 14,271 14,271 17,267 17,267 18,266 18,266 20,263 20,263 18,262 18,262 17,259 17,259 16,253 16,253 17,252 17,252 18 [...]
+MULTIPOLYGON (((323 0,367 0,367 2,365 2,365 3,363 3,363 4,362 4,362 5,363 5,363 7,365 7,365 6,368 6,368 5,370 5,370 4,373 4,373 3,378 3,378 4,380 4,380 5,381 5,381 6,382 6,382 7,383 7,383 8,385 8,385 9,387 9,387 10,390 10,390 9,395 9,395 8,393 8,393 7,392 7,392 3,393 3,393 2,392 2,392 1,393 1,393 0,503 0,503 12,502 12,502 13,501 13,501 14,499 14,499 15,498 15,498 16,497 16,497 17,496 17,496 18,493 18,493 19,492 19,492 20,498 20,498 19,502 19,502 18,503 18,503 33,501 33,501 32,491 32,491 [...]
+MULTIPOLYGON (((426 13,429 13,429 14,430 14,430 13,431 13,431 14,438 14,438 15,441 15,441 14,443 14,443 15,444 15,444 16,446 16,446 17,449 17,449 18,455 18,455 19,457 19,457 20,459 20,459 21,460 21,460 20,463 20,463 21,465 21,465 23,466 23,466 27,467 27,467 29,466 29,466 30,465 30,465 32,466 32,466 33,472 33,472 34,473 34,473 35,479 35,479 36,481 36,481 35,485 35,485 34,489 34,489 33,491 33,491 32,501 32,501 33,503 33,503 55,501 55,501 56,500 56,500 57,493 57,493 58,491 58,491 61,492 61, [...]
diff --git a/src/tests/good_test1_features.dbf b/src/tests/good_test1_features.dbf
new file mode 100644
index 0000000..57a67c8
Binary files /dev/null and b/src/tests/good_test1_features.dbf differ
diff --git a/src/tests/good_test1_features.shp b/src/tests/good_test1_features.shp
new file mode 100644
index 0000000..3b469b9
Binary files /dev/null and b/src/tests/good_test1_features.shp differ
diff --git a/src/tests/good_test1_features.shx b/src/tests/good_test1_features.shx
new file mode 100644
index 0000000..7a0f291
Binary files /dev/null and b/src/tests/good_test1_features.shx differ
diff --git a/src/tests/good_test1_features.wkt b/src/tests/good_test1_features.wkt
new file mode 100644
index 0000000..4dac323
--- /dev/null
+++ b/src/tests/good_test1_features.wkt
@@ -0,0 +1,20 @@
+MULTIPOLYGON (((411 204,412 204,412 205,411 205,411 204)),((413 205,414 205,414.0 205.9,413.9 206.0,413 206,413 205)),((258 206,259 206,259 208,258.1 208.0,258.0 207.9,258 206)),((414 206,417 206,417 207,414 207,414 206)),((257 208,258 208,258 209,257 209,257 208)),((466 248,467 248,467 252,466 252,466 248)),((43 252,44 252,44 254,43 254,43 252)),((66 260,67 260,67 262,66 262,66 260)),((362 355,363 355,363 356,362 356,362 355)),((465 359,466 359,466 360,465 360,465 359)),((485 381,486 38 [...]
+MULTIPOLYGON (((37 0,42 0,42 2,43 2,43.0 3.9,42.9 4.0,41 4,41 3,40 3,40 2,38 2,38 1,37 1,37 0)),((57 0,63 0,63 2,64 2,64 4,65 4,65 5,66 5,66 7,67 7,67 8,68 8,68 9,64 9,64 8,61 8,61 6,60 6,60 5,59 5,59 3,58 3,58 1,57 1,57 0)),((114 0,116 0,116 3,117 3,117 7,118 7,118 8,117 8,117 12,116 12,116 15,115 15,115 17,114 17,114 19,113 19,113 20,112 20,112 21,111 21,111 22,110 22,110 24,109 24,109 25,108 25,108 26,106 26,106 27,104 27,104 28,98 28,98 26,100 26,100 25,102 25,102 24,103 24,103 23,10 [...]
+MULTIPOLYGON (((50 0,54 0,54 1,55 1,55 5,56 5,56 6,57 6,57 7,55 7,55 6,54 6,54 5,53 5,53 4,52 4,52 3,51 3,51 1,50 1,50 0)),((77 0,78 0,78.0 0.9,77.9 1.0,77 1,77 0)),((99 0,105 0,105 2,104 2,104 3,98 3,98 1,99 1,99 0)),((209 0,212 0,212 3,211 3,211 2,209 2,209 0)),((256 0,258 0,258 1,259 1,259 3,257 3,257 2,256 2,256 0)),((357 0,361 0,361 1,362 1,362 3,361 3,361 4,360 4,360 3,359 3,359 2,358 2,358 1,357 1,357 0)),((78 1,79 1,79 2,78 2,78 1)),((121 10,122 10,122 11,121.1 11.0,121.0 10.9,12 [...]
+MULTIPOLYGON (((31 0,36 0,36 1,35 1,35 2,32 2,32 1,31 1,31 0)),((54 0,57 0,57 1,58 1,58 3,59 3,59 5,60 5,60 6,61 6,61 8,62 8,62 9,59 9,59 8,57 8,57 6,56 6,56 5,55 5,55 1,54 1,54 0)),((78 0,82 0,82 2,81 2,81 3,79 3,79 1,78 1,78 0)),((95 0,98 0,98 4,97 4,97 5,94 5,94 6,93 6,93 7,89 7,89 5,90 5,90 4,91 4,91 3,93 3,93 2,94 2,94 1,95 1,95 0)),((139 0,144 0,144 1,143 1,143 6,142 6,142 7,141 7,141 9,139 9,139 10,137 10,137 11,135 11,135 10,136 10,136 9,137 9,137 7,138 7,138 4,139 4,139 0)),((21 [...]
+MULTIPOLYGON (((24 0,29 0,29 1,30 1,30 2,28 2,28 3,26 3,26 2,25 2,25 1,24 1,24 0)),((72 0,77 0,77 1,78 1,78 2,79 2,79 3,80 3,80 4,81 4,81 5,82 5,82.0 5.9,81.9 6.0,81 6,81.0 6.9,80.9 7.0,79 7,79 6,78 6,78 5,76 5,76 4,74 4,74 2,73 2,73 1,72 1,72 0)),((98 0,99 0,99 1,98 1,98 0)),((124 0,126 0,126 4,125 4,125 7,124 7,124 12,123 12,123 15,122 15,122 16,121 16,121 18,120 18,120 20,119 20,119 21,118 21,118 22,116 22,116 21,117 21,117 19,118 19,118 17,119 17,119 15,120 15,120 13,121 13,121 11,12 [...]
+MULTIPOLYGON (((0 0,5 0,5 1,6 1,6 0,9 0,9 1,13 1,13 2,15 2,15 3,16 3,16 4,18 4,18 5,19 5,19 6,17 6,17 7,14 7,14 8,7 8,7 7,6 7,6 6,5 6,5 5,3 5,3 4,0 4,0 0)),((47 0,50 0,50 1,51 1,51 3,52 3,52 5,53 5,53 7,54 7,54 8,52 8,52 7,51 7,51 6,50 6,50 4,49 4,49 3,48 3,48 2,47 2,47 0)),((66 0,72 0,72 1,73 1,73 2,74 2,74 4,69 4,69 3,68 3,68 2,66 2,66 0)),((105 0,109 0,109 2,110 2,110 3,109 3,109 4,108 4,108 5,107 5,107 6,105 6,105 7,103 7,103 4,104 4,104 2,105 2,105 0)),((126 0,131 0,131 1,130 1,130 [...]
+MULTIPOLYGON (((36 0,37 0,37 1,38 1,38 2,40 2,40 3,41 3,41 4,43 4,43 5,44 5,44 7,45 7,45 9,46 9,46 11,47 11,47 12,48 12,48 14,49 14,49 16,50 16,50 17,51 17,51 19,52 19,52 20,53 20,53 22,54 22,54 23,55 23,55 25,56 25,56 26,57 26,57 27,58 27,58 28,59 28,59 30,60 30,60 31,61 31,61 32,62 32,62 33,63 33,63 34,64 34,64 35,65 35,65 36,66 36,66 37,67 37,67 38,68 38,68 39,69 39,69 40,70 40,70 41,72 41,72 42,73 42,73 43,75 43,75 45,76 45,76 51,75 51,75 50,74 50,74 49,72 49,72 48,71 48,71 47,69 47, [...]
+MULTIPOLYGON (((82 0,95 0,95 1,94 1,94 2,93 2,93 3,91 3,91 4,90 4,90 5,89 5,89 6,87 6,87 5,85 5,85 4,83 4,83 2,82 2,82 0)),((119 0,122 0,122 3,121 3,121 5,119 5,119 4,118 4,118 3,119 3,119 0)),((195 0,198 0,198 3,196 3,196 4,195.1 4.0,195.0 3.9,195 0)),((275 0,280 0,280 1,281 1,281 2,282 2,282 3,280 3,280 4,278 4,278 3,277 3,277 2,276 2,276 1,275 1,275 0)),((320 0,322 0,322 1,321 1,321 2,320 2,320 0)),((420 0,424 0,424 1,425 1,425 2,426 2,426 3,427 3,427 4,428 4,428 5,429 5,429 6,430 6,4 [...]
+MULTIPOLYGON (((29 0,31 0,31 1,32 1,32 9,34 9,34 10,33 10,33 11,31 11,31 12,30 12,30.0 13.9,29.9 14.0,29 14,29 15,27 15,27 16,24 16,24 15,22 15,22 14,21 14,21 11,23 11,23 10,25 10,25 9,26.9 9.0,27.0 9.1,27 11,28 11,28 10,29 10,29 9,27 9,27 7,28 7,28 6,27 6,27 5,28 5,28 2,30 2,30 1,29 1,29 0)),((122 0,124 0,124 2,123 2,123 4,122 4,122 7,121 7,121 3,122 3,122 0)),((144 0,145 0,145 1,144.1 1.0,144.0 0.9,144 0)),((193 0,195 0,195 4,194 4,194 15,193 15,193 18,192 18,192 21,190 21,190 22,189.1 [...]
+MULTIPOLYGON (((5 0,6 0,6 1,5 1,5 0)),((9 0,11 0,11 1,9 1,9 0)),((42 0,47 0,47 3,48 3,48 4,47 4,47.0 4.9,46.9 5.0,46 5,46 4,45 4,45 3,44 3,44 2,42 2,42 0)),((63 0,66 0,66 2,68 2,68 3,69 3,69 4,70 4,70 5,71 5,71 6,70 6,70 7,66 7,66 5,65 5,65 4,64 4,64 2,63 2,63 0)),((109 0,113 0,113 1,114 1,114 2,113 2,113 7,112 7,112 10,111 10,111 13,110 13,110 15,108 15,108 14,107 14,107 11,108 11,108 10,107 10,107 8,108 8,108 4,109 4,109 3,110 3,110 2,109 2,109 0)),((131 0,134 0,134 3,133 3,133 6,132 6 [...]
+MULTIPOLYGON (((23 0,24 0,24 1,25 1,25 2,26 2,26 3,28 3,28 5,27 5,27 6,28 6,28 7,27 7,27.0 8.9,26.9 9.0,25 9,25 10,23 10,23 11,22 11,22 10,21 10,21 9,20 9,20 7,19 7,19 5,21 5,21 3,22 3,22 2,23 2,23 0)),((231 0,239 0,239 1,238 1,238 2,237 2,237 3,236 3,236 4,237 4,237.0 8.9,236.9 9.0,236 9,236 11,237 11,237 9,238 9,238 13,237 13,237 14,239 14,239 18,240 18,240 21,241 21,241 26,240 26,240 24,239 24,239 23,238 23,238 22,237 22,237 20,236 20,236 16,235 16,235 14,234 14,234 12,233 12,233 10,2 [...]
+MULTIPOLYGON (((11 0,23 0,23 2,22 2,22 3,21 3,21 5,18 5,18 4,16 4,16 3,15 3,15 2,13 2,13 1,11 1,11 0)),((182 0,187 0,187 1,188 1,188 2,189 2,189 5,184 5,184 2,182 2,182 0)),((226 0,231 0,231 2,230 2,230 3,231 3,231 8,230 8,230 7,226 7,226 0)),((287 0,293 0,293 1,294 1,294 6,293 6,293 7,294 7,294 8,295 8,295 10,294 10,294 11,293 11,293 12,291 12,291 13,290 13,290 9,291 9,291 8,292 8,292 7,289 7,289 4,288 4,288 1,287 1,287 0)),((349 0,352 0,352 1,353 1,353 2,354 2,354 4,355 4,355 7,356 7,3 [...]
+MULTIPOLYGON (((190 0,193 0,193 1,190 1,190 0)),((327 0,333 0,333 2,329 2,329 1,327 1,327 0)),((425 0,429 0,429 3,428 3,428 2,426 2,426 1,425 1,425 0)),((238 1,239 1,239 2,238.1 2.0,238.0 1.9,238 1)),((388 1,389 1,389 2,388.1 2.0,388.0 1.9,388 1)),((237 2,238 2,238.0 3.9,237.9 4.0,236 4,236 3,237 3,237 2)),((387 2,388 2,388 4,387 4,387 6,386 6,386 7,385 7,385 8,384 8,384 10,383 10,383 11,382 11,382 12,381 12,381 13,378 13,378 12,379 12,379 10,380 10,380 9,382 9,382 8,383 8,383 7,384 7,38 [...]
+MULTIPOLYGON (((390 0,403 0,403 1,402 1,402 2,400 2,400 3,398 3,398 4,396 4,396 5,395 5,395 6,394 6,394 7,393 7,393 8,391 8,391 9,387 9,387 8,386 8,386 7,387 7,387 4,388 4,388 2,389 2,389 1,390 1,390 0)),((39 8,40 8,40 9,41 9,41 10,42 10,42 11,43 11,43 13,44 13,44 14,45 14,45 15,46 15,46 16,47 16,47 18,48 18,48 21,49 21,49 22,48 22,48 24,46 24,46 23,42 23,42 22,41 22,41 21,40 21,40 20,37 20,37 13,36 13,36 12,37 12,37 11,38 11,38.0 10.1,38.1 10.0,39 10,39 8)),((490 8,500 8,500 13,494 13,4 [...]
+MULTIPOLYGON (((187 0,190 0,190 1,189 1,189 2,188 2,188 1,187 1,187 0)),((338 0,349 0,349 3,350 3,350 4,351 4,351 5,352 5,352 6,353 6,353 7,354 7,354 8,355 8,355 10,356 10,356 12,357 12,357 13,356 13,356 14,357 14,357 18,356 18,356 22,355 22,355 24,354 24,354 27,353 27,353 29,352 29,352 31,351 31,351 32,350 32,350 33,347 33,347 32,346 32,346 30,345 30,345 29,346 29,346 27,347 27,347 24,348 24,348 22,349 22,349 13,348 13,348 10,347 10,347 9,345 9,345 8,344 8,344 7,343 7,343 6,342 6,342 5, [...]
+MULTIPOLYGON (((113 0,114 0,114 1,113 1,113 0)),((163 0,180 0,180 3,179 3,179 6,178 6,178 8,177 8,177 9,175 9,175 10,172 10,172 9,171 9,171 8,169 8,169 6,168 6,168 5,167 5,167 4,166 4,166 3,165 3,165 2,164 2,164 1,163 1,163 0)),((268 0,269 0,269.0 0.9,268.9 1.0,268 1,268 0)),((269 1,270 1,270.0 1.9,269.9 2.0,269 2,269 1)),((405 1,409 1,409 6,408 6,408 7,407 7,407 6,406 6,406 7,405 7,405 6,404.1 6.0,404.0 5.9,404 3,405 3,405 1)),((475 1,481 1,481 2,485 2,485 3,498 3,498 2,500 2,500 6,486 [...]
+MULTIPOLYGON (((385 0,390 0,390 1,388 1,388 2,387 2,387 3,386 3,386 5,385 5,385 6,384 6,384 7,383 7,383 8,382 8,382 9,380.1 9.0,380.0 8.9,380 6,381 6,381 5,382 5,382 4,383 4,383 2,384 2,384 1,385 1,385 0)),((379 9,380 9,380 10,379 10,379 9)),((0 10,2 10,2 11,3 11,3 12,4 12,4 16,3 16,3 17,0 17,0 10)),((211 11,212 11,212 12,213 12,213 13,214 13,214 16,215 16,215 22,214 22,214 20,213 20,213 18,212 18,212 15,211 15,211 11)),((247 13,248 13,248 14,250 14,250 15,251 15,251 17,250 17,250 20,246 [...]
+MULTIPOLYGON (((382 0,385 0,385 1,384 1,384 2,383 2,383 4,382 4,382 5,381 5,381 6,379 6,379 5,380 5,380 3,381 3,381 1,382 1,382 0)),((495 13,500 13,500 20,499 20,499 21,493 21,493 23,491 23,491 24,490 24,490 23,489.1 23.0,489.0 22.9,489 21,490 21,490 17,491 17,491 16,492 16,492 15,493 15,493 14,495 14,495 13)),((211 15,212 15,212.0 17.9,211.9 18.0,211 18,211 15)),((61 18,62 18,62 19,65 19,65 20,66 20,66 21,67 21,67 22,68 22,68 23,69 23,69 25,70 25,70 28,71 28,71 30,70 30,70 29,68 29,68 2 [...]
+MULTIPOLYGON (((377 0,382 0,382 1,381 1,381 3,380 3,380 5,379 5,379 6,378 6,378 7,375 7,375 6,376 6,376 3,377 3,377 0)),((401 7,403 7,403 8,402 8,402 10,400 10,400 11,399 11,399 12,396 12,396 11,397 11,397 10,398 10,398 9,399 9,399 8,401 8,401 7)),((207 12,208 12,208 13,210 13,210 14,211 14,211 18,212 18,212 23,211 23,211 26,212 26,212 27,213 27,213 29,214 29,214 34,213 34,213 32,212 32,212 31,211 31,211 29,210 29,210 26,209 26,209 24,208 24,208 19,207 19,207 12)),((422 15,424 15,424 16, [...]
+MULTIPOLYGON (((373 0,377 0,377 3,376 3,376 6,373 6,373 5,372 5,372 3,373 3,373 0)),((402 1,405 1,405 3,404 3,404 6,403 6,403 7,401 7,401 8,399 8,399 9,398 9,398 10,397 10,397 11,396 11,396 12,391 12,391 10,392 10,392 9,393 9,393 7,394 7,394 6,395 6,395 5,396 5,396 4,398 4,398 3,400 3,400 2,402 2,402 1)),((204 12,207 12,207 19,208 19,208 24,209 24,209 26,210 26,210 29,211 29,211 31,212 31,212 32,213 32,213 36,212 36,212 37,210 37,210 38,208 38,208 39,207 39,207 38,206 38,206 35,205 35,20 [...]
diff --git a/src/tests/good_test1_gradient_ndv.pbm b/src/tests/good_test1_gradient_ndv.pbm
new file mode 100644
index 0000000..c20f132
Binary files /dev/null and b/src/tests/good_test1_gradient_ndv.pbm differ
diff --git a/src/tests/good_test1_gradient_ndv_inv.pbm b/src/tests/good_test1_gradient_ndv_inv.pbm
new file mode 100644
index 0000000..cf4c4b6
Binary files /dev/null and b/src/tests/good_test1_gradient_ndv_inv.pbm differ
diff --git a/src/tests/good_test1_gradient_valid.pbm b/src/tests/good_test1_gradient_valid.pbm
new file mode 100644
index 0000000..cf4c4b6
Binary files /dev/null and b/src/tests/good_test1_gradient_valid.pbm differ
diff --git a/src/tests/good_test1_maze.ppm b/src/tests/good_test1_maze.ppm
new file mode 100644
index 0000000..ce5b8ba
Binary files /dev/null and b/src/tests/good_test1_maze.ppm differ
diff --git a/src/tests/good_test1_nan1.pbm b/src/tests/good_test1_nan1.pbm
new file mode 100644
index 0000000..921bdcb
Binary files /dev/null and b/src/tests/good_test1_nan1.pbm differ
diff --git a/src/tests/good_test1_nan2.pbm b/src/tests/good_test1_nan2.pbm
new file mode 100644
index 0000000..96650ae
Binary files /dev/null and b/src/tests/good_test1_nan2.pbm differ
diff --git a/src/tests/good_test1_nan3.pbm b/src/tests/good_test1_nan3.pbm
new file mode 100644
index 0000000..921bdcb
Binary files /dev/null and b/src/tests/good_test1_nan3.pbm differ
diff --git a/src/tests/good_test1_nan4.pbm b/src/tests/good_test1_nan4.pbm
new file mode 100644
index 0000000..28803a7
Binary files /dev/null and b/src/tests/good_test1_nan4.pbm differ
diff --git a/src/tests/good_test1_noise.ppm b/src/tests/good_test1_noise.ppm
new file mode 100644
index 0000000..a4aee7f
Binary files /dev/null and b/src/tests/good_test1_noise.ppm differ
diff --git a/src/tests/good_test1_noise_dp3.ppm b/src/tests/good_test1_noise_dp3.ppm
new file mode 100644
index 0000000..150b599
Binary files /dev/null and b/src/tests/good_test1_noise_dp3.ppm differ
diff --git a/src/tests/good_test2_1.ppm b/src/tests/good_test2_1.ppm
index 2bacda6..d99bda4 100644
Binary files a/src/tests/good_test2_1.ppm and b/src/tests/good_test2_1.ppm differ
diff --git a/src/tests/good_test2_2.ppm b/src/tests/good_test2_2.ppm
index e308f38..d0263e9 100644
Binary files a/src/tests/good_test2_2.ppm and b/src/tests/good_test2_2.ppm differ
diff --git a/src/tests/good_test2_3.ppm b/src/tests/good_test2_3.ppm
index e5be983..91f213e 100644
Binary files a/src/tests/good_test2_3.ppm and b/src/tests/good_test2_3.ppm differ
diff --git a/src/tests/good_test2_4.ppm b/src/tests/good_test2_4.ppm
index a5e04c1..c0be0d7 100644
Binary files a/src/tests/good_test2_4.ppm and b/src/tests/good_test2_4.ppm differ
diff --git a/src/tests/good_test2_5.ppm b/src/tests/good_test2_5.ppm
index 09b8d08..78a6bef 100644
Binary files a/src/tests/good_test2_5.ppm and b/src/tests/good_test2_5.ppm differ
diff --git a/src/tests/good_test2_maze.ppm b/src/tests/good_test2_maze.ppm
index 935186a..ce5b8ba 100644
Binary files a/src/tests/good_test2_maze.ppm and b/src/tests/good_test2_maze.ppm differ
diff --git a/src/tests/good_test2_noise.ppm b/src/tests/good_test2_noise.ppm
index a133861..a4aee7f 100644
Binary files a/src/tests/good_test2_noise.ppm and b/src/tests/good_test2_noise.ppm differ
diff --git a/src/tests/good_test2_noise_dp3.ppm b/src/tests/good_test2_noise_dp3.ppm
index 59bd360..150b599 100644
Binary files a/src/tests/good_test2_noise_dp3.ppm and b/src/tests/good_test2_noise_dp3.ppm differ
diff --git a/src/tests/good_test3_dem.tif b/src/tests/good_test3_dem.tif
new file mode 100644
index 0000000..d18a946
Binary files /dev/null and b/src/tests/good_test3_dem.tif differ
diff --git a/src/tests/good_test2_4.ppm b/src/tests/good_test3_histeq.tif
similarity index 53%
copy from src/tests/good_test2_4.ppm
copy to src/tests/good_test3_histeq.tif
index a5e04c1..cb1761e 100644
Binary files a/src/tests/good_test2_4.ppm and b/src/tests/good_test3_histeq.tif differ
diff --git a/src/tests/test1.sh b/src/tests/test1.sh
index 1adda41..6ca6947 100755
--- a/src/tests/test1.sh
+++ b/src/tests/test1.sh
@@ -1,4 +1,4 @@
-#!/bin/sh
+#!/bin/bash
rm -f out_test1_*
@@ -15,8 +15,11 @@ $BINDIR/gdal_trace_outline testcase_maze.png -ndv 255 -out-cs xy -wkt-out out_t
$BINDIR/gdal_trace_outline testcase_noise.png -b 1 -ndv 0 -out-cs xy -wkt-out out_test1_noise.wkt -report out_test1_noise.ppm -split-polys -dp-toler 0
$BINDIR/gdal_trace_outline testcase_noise.png -b 1 -ndv 0 -out-cs xy -wkt-out out_test1_noise_dp3.wkt -report out_test1_noise_dp3.ppm -split-polys -dp-toler 3
-$BINDIR/gdal_trace_outline testcase_3.tif -out-cs xy -wkt-out out_test1_3_classify.wkt -dp-toler 0 -classify
-$BINDIR/gdal_trace_outline pal.tif -out-cs xy -wkt-out out_test1_3_classify_pal.wkt -dp-toler 0 -classify
+$BINDIR/gdal_trace_outline testcase_3.tif -out-cs xy -wkt-out out_test1_3_classify.wkt -ogr-out out_test1_3_classify.shp -dp-toler 0 -classify
+$BINDIR/gdal_trace_outline testcase_3_paletted.tif -out-cs xy -wkt-out out_test1_3_classify_pal.wkt -ogr-out out_test1_3_classify_pal.shp -dp-toler 0 -classify
+$BINDIR/gdal_trace_outline testcase_features.png -out-cs xy -wkt-out out_test1_features.wkt -ogr-out out_test1_features.shp -dp-toler 0 -classify
+$BINDIR/gdal_trace_outline testcase_double.tif -out-cs xy -wkt-out out_test1_double.wkt -ogr-out out_test1_double.shp -dp-toler 0 -classify
+$BINDIR/gdal_trace_outline testcase_double.tif -out-cs xy -wkt-out out_test1_double_clip.wkt -dp-toler 0 -classify -valid-range '3..6'
$BINDIR/gdal_list_corners -inspect-rect4 -erosion -ndv 0 testcase_4.png -report out_test1_4-rect.ppm > out_test1_4-rect.wkt
@@ -27,10 +30,93 @@ $BINDIR/gdal_get_projected_bounds -s_wkt good_test1_1_en.wkt -s_srs '+proj=utm +
$BINDIR/gdal_list_corners testcase_1.tif > out_test1_1_metdata.yml
$BINDIR/gdal_make_ndv_mask -ndv '155 52 52' -ndv '24 173 79' testcase_3.tif out_test1_3_ndvmask.pbm
-$BINDIR/gdal_make_ndv_mask -ndv '155 52 52' -ndv '24 173 79..81' testcase_3.tif out_test1_3_ndvmask2.pbm
+$BINDIR/gdal_make_ndv_mask -ndv '155 52 52' -ndv '24 173 79.9..80.1' testcase_3.tif out_test1_3_ndvmask2.pbm
+
+# Make a gradient image.
+python <<END
+import numpy as np
+import osgeo.gdal as gdal
+
+driver = gdal.GetDriverByName('GTiff')
+
+(w, h) = (256, 256)
+
+x = np.outer(np.ones(h), range(w))
+y = np.outer(range(h), np.ones(w))
+
+dst_ds = driver.Create('gradient1.tif', w, h, 2, gdal.GDT_Byte)
+dst_ds.GetRasterBand(1).WriteArray(x, 0, 0)
+dst_ds.GetRasterBand(2).WriteArray(y, 0, 0)
+dst_ds = None
+
+dst_ds = driver.Create('gradient2.tif', w, h, 1, gdal.GDT_Float64)
+dst_ds.GetRasterBand(1).WriteArray((y+1)/(x+1), 0, 0)
+dst_ds = None
+END
+gdal_merge_vrt -in gradient1.tif -in gradient2.tif -out gradient3.tif
+
+$BINDIR/gdal_make_ndv_mask \
+ -ndv '10..30 30..70 *' \
+ -ndv '* * 4..Inf' \
+ -ndv '100..140 50..80 0.5..Inf' \
+ -ndv '* * 0.8..0.3' \
+ -ndv '* * 0.3..0.4' \
+ gradient3.tif out_test1_gradient_ndv.pbm
+
+$BINDIR/gdal_make_ndv_mask \
+ -ndv '10..30 30..70 *' \
+ -ndv '* * 4..Inf' \
+ -ndv '100..140 50..80 0.5..Inf' \
+ -ndv '* * 0.8..0.3' \
+ -ndv '* * 0.3..0.4' \
+ -invert \
+ gradient3.tif out_test1_gradient_ndv_inv.pbm
+
+$BINDIR/gdal_make_ndv_mask \
+ -valid-range '10..30 30..70 *' \
+ -valid-range '* * 4..Inf' \
+ -valid-range '100..140 50..80 0.5..Inf' \
+ -valid-range '* * 0.8..0.3' \
+ -valid-range '* * 0.3..0.4' \
+ gradient3.tif out_test1_gradient_valid.pbm
+
+python <<END
+import numpy as np
+import osgeo.gdal as gdal
+
+driver = gdal.GetDriverByName('GTiff')
+(w, h) = (256, 256)
+dst_ds = driver.Create('has_nan.tif', w, h, 2, gdal.GDT_Float64)
+
+vals = np.zeros((w, h))
+vals[30:40,100:110] = np.nan
+vals[60:70,100:110] = np.inf
+dst_ds.GetRasterBand(1).WriteArray(vals, 0, 0)
+
+vals = np.zeros((w, h))
+vals[30:40,150:160] = np.nan
+vals[60:70,150:160] = -np.inf
+dst_ds.GetRasterBand(2).WriteArray(vals, 0, 0)
+
+dst_ds = None
+END
+
+# NDV is specified here just to make the script run. But the point is to see that NaN values
+# (in either band) are considered to be NDV.
+$BINDIR/gdal_make_ndv_mask has_nan.tif out_test1_nan1.pbm -ndv 7
+# Check detection of infinity
+$BINDIR/gdal_make_ndv_mask has_nan.tif out_test1_nan2.pbm -ndv 'Inf *'
+$BINDIR/gdal_make_ndv_mask has_nan.tif out_test1_nan3.pbm -ndv '* Inf'
+$BINDIR/gdal_make_ndv_mask has_nan.tif out_test1_nan4.pbm -ndv '* -Inf'
echo '####################'
+for i in out_test1_* ; do
+ if [ ! -e ${i/out/good} ] ; then
+ echo "!!! ${i/out/good} doesn't exist"
+ fi
+done
+
for i in good_test1_* ; do
if diff --brief $i ${i/good/out} ; then
echo "GOOD ${i/good_/}"
diff --git a/src/tests/test2.sh b/src/tests/test2.sh
index aa163c4..1111277 100755
--- a/src/tests/test2.sh
+++ b/src/tests/test2.sh
@@ -1,4 +1,4 @@
-#!/bin/sh
+#!/bin/bash
rm -f out_test2_*
@@ -16,6 +16,12 @@ $BINDIR/gdal_trace_outline testcase_noise.png -b 1 -ndv 0 -out-cs xy -wkt-out
echo '####################'
+for i in out_test2_* ; do
+ if [ ! -e ${i/out/good} ] ; then
+ echo "!!! ${i/out/good} doesn't exist"
+ fi
+done
+
for i in good_test2_* ; do
if diff --brief $i ${i/good/out} ; then
echo "GOOD ${i/good_/}"
diff --git a/src/tests/test3.sh b/src/tests/test3.sh
index c098edc..b748c12 100755
--- a/src/tests/test3.sh
+++ b/src/tests/test3.sh
@@ -1,4 +1,4 @@
-#!/bin/sh
+#!/bin/bash
rm -f out_test3_test3_*
@@ -10,6 +10,12 @@ $BINDIR/gdal_contrast_stretch -ndv '0..255 0..255 0..255 0' -histeq 50 testcase_
echo '####################'
+for i in out_test3_* ; do
+ if [ ! -e ${i/out/good} ] ; then
+ echo "!!! ${i/out/good} doesn't exist"
+ fi
+done
+
for i in good_test3_* ; do
if diff --brief $i ${i/good/out} ; then
echo "GOOD ${i/good_/}"
diff --git a/src/tests/pal.tif b/src/tests/testcase_3_paletted.tif
similarity index 100%
rename from src/tests/pal.tif
rename to src/tests/testcase_3_paletted.tif
diff --git a/src/tests/testcase_double.tif b/src/tests/testcase_double.tif
new file mode 100644
index 0000000..d9de57a
Binary files /dev/null and b/src/tests/testcase_double.tif differ
diff --git a/src/tests/testcase_features.png b/src/tests/testcase_features.png
new file mode 100644
index 0000000..a576cdf
Binary files /dev/null and b/src/tests/testcase_features.png differ
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/pkg-grass/dans-gdal-scripts.git
More information about the Pkg-grass-devel
mailing list