[geos-commits] [SCM] GEOS branch master updated. 352d07a45ea4d7b8a17962e7d02ba9a69dce891d
git at osgeo.org
git at osgeo.org
Thu Nov 26 10:18:58 PST 2020
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "GEOS".
The branch, master has been updated
via 352d07a45ea4d7b8a17962e7d02ba9a69dce891d (commit)
from 34442d22bacbdc733d41d6948f8f3e33044e68dc (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
commit 352d07a45ea4d7b8a17962e7d02ba9a69dce891d
Author: Paul Ramsey <pramsey at cleverelephant.ca>
Date: Thu Nov 26 10:18:02 2020 -0800
Add MortonCode and HilbertCode for future use.
Also add relevant tests. Back out removal of
the float-store flag for 32-bit platforms.
diff --git a/configure.ac b/configure.ac
index fceccc9..9b64ba3 100644
--- a/configure.ac
+++ b/configure.ac
@@ -202,8 +202,17 @@ AC_LIBTOOL_COMPILER_OPTION([if $compiler supports -pedantic], [dummy_cv_pedantic
AC_LIBTOOL_COMPILER_OPTION([if $compiler supports -Wall], [dummy_cv_wall], [-Wall], [], [WARNFLAGS="$WARNFLAGS -Wall"], [])
AC_LIBTOOL_COMPILER_OPTION([if $compiler supports -Wno-long-long], [dummy_cv_wno_long_long], [-Wno-long-long], [], [WARNFLAGS="$WARNFLAGS -Wno-long-long"], [])
+# To make numerical computation more stable, we use --ffloat-store
+# on 32-bit platforms
+NUMERICFLAGS=""
+# 32-bit platforms have a 4-byte pointer
+AC_CHECK_SIZEOF([void *])
+if test $ac_cv_sizeof_void_p -eq 4 ; then
+ AC_LIBTOOL_COMPILER_OPTION([if $compiler supports -ffloat-store], [dummy_cv_ffloat_store], [-ffloat-store], [], [NUMERICFLAGS="$NUMERICFLAGS -ffloat-store"], [])
+fi
+
HUSHWARNING="-DUSE_UNSTABLE_GEOS_CPP_API"
-DEFAULTFLAGS="${WARNFLAGS} ${HUSHWARNING} ${OVERLAYNG_FLAGS}"
+DEFAULTFLAGS="${WARNFLAGS} ${NUMERICFLAGS} ${HUSHWARNING} ${OVERLAYNG_FLAGS}"
AM_CXXFLAGS="${AM_CXXFLAGS} ${DEFAULTFLAGS}"
AM_CFLAGS="${AM_CFLAGS} ${DEFAULTFLAGS}"
@@ -442,6 +451,8 @@ AC_OUTPUT([
include/geos/simplify/Makefile
include/geos/triangulate/Makefile
include/geos/triangulate/quadedge/Makefile
+ include/geos/shape/Makefile
+ include/geos/shape/fractal/Makefile
include/geos/util/Makefile
include/geos/version.h
src/index/Makefile
@@ -476,6 +487,8 @@ AC_OUTPUT([
src/simplify/Makefile
src/triangulate/Makefile
src/triangulate/quadedge/Makefile
+ src/shape/Makefile
+ src/shape/fractal/Makefile
src/util/Makefile
swig/geos.i
swig/Makefile
diff --git a/include/geos/Makefile.am b/include/geos/Makefile.am
index 1b2094c..09ba8db 100644
--- a/include/geos/Makefile.am
+++ b/include/geos/Makefile.am
@@ -16,6 +16,7 @@ SUBDIRS = \
precision \
simplify \
triangulate \
+ shape \
util
EXTRA_DIST =
diff --git a/include/geos/shape/Makefile.am b/include/geos/shape/Makefile.am
new file mode 100644
index 0000000..1c303f2
--- /dev/null
+++ b/include/geos/shape/Makefile.am
@@ -0,0 +1,11 @@
+#
+# This file is part of project GEOS (http://trac.osgeo.org/geos/)
+#
+SUBDIRS = \
+ fractal
+
+EXTRA_DIST =
+
+geosdir = $(includedir)/geos/shape
+
+geos_HEADERS =
diff --git a/include/geos/shape/fractal/HilbertCode.h b/include/geos/shape/fractal/HilbertCode.h
new file mode 100644
index 0000000..46012ef
--- /dev/null
+++ b/include/geos/shape/fractal/HilbertCode.h
@@ -0,0 +1,127 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2020 Paul Ramsey <pramsey at cleverelephant.ca>
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU Lesser General Public Licence as published
+ * by the Free Software Foundation.
+ * See the COPYING file for more information.
+ *
+ **********************************************************************/
+
+
+#pragma once
+
+#include <geos/export.h>
+#include <string>
+
+// Forward declarations
+namespace geos {
+namespace geom {
+class Coordinate;
+}
+}
+
+namespace geos {
+namespace shape { // geos.shape
+namespace fractal { // geos.shape.fractal
+
+/**
+ * Encodes points as the index along finite planar Hilbert curves.
+ *
+ * The planar Hilbert Curve is a continuous space-filling curve.
+ * In the limit the Hilbert curve has infinitely many vertices and fills
+ * the space of the unit square.
+ * A sequence of finite approximations to the infinite Hilbert curve
+ * is defined by the level number.
+ * The finite Hilbert curve at level n H(n) contains 2^n+1 points.
+ * Each finite Hilbert curve defines an ordering of the
+ * points in the 2-dimensional range square containing the curve.
+ * Curves fills the range square of side 2^level.
+ * Curve points have ordinates in the range [0, 2^level - 1].
+ * The index of a point along a Hilbert curve is called the Hilbert code.
+ * The code for a given point is specific to the level chosen.
+ *
+ * This implementation represents codes using 32-bit integers.
+ * This allows levels 0 to 16 to be handled.
+ * The class supports encoding points in the range of a given level curve
+ * and decoding the point for a given code value.
+ *
+ * The Hilbert order has the property that it tends to preserve locality.
+ * This means that codes which are near in value will have spatially proximate
+ * points. The converse is not always true - the delta between
+ * codes for nearby points is not always small. But the average delta
+ * is small enough that the Hilbert order is an effective way of linearizing space
+ * to support range queries.
+ *
+ * @author Martin Davis
+ *
+ * @see MortonCode
+ */
+class GEOS_DLL HilbertCode {
+
+public:
+
+ /**
+ * The maximum curve level that can be represented.
+ */
+ static constexpr int MAX_LEVEL = 16;
+
+ static geom::Coordinate decode(uint32_t level, uint32_t i);
+
+ static uint32_t encode(uint32_t level, uint32_t x, uint32_t y);
+
+ /**
+ * The number of points in the curve for the given level.
+ * The number of points is 2^(2 * level).
+ *
+ * @param level the level of the curve
+ * @return the number of points
+ */
+ static uint32_t levelSize(uint32_t level);
+
+ /**
+ * The maximum ordinate value for points
+ * in the curve for the given level.
+ * The maximum ordinate is 2^level - 1.
+ *
+ * @param level the level of the curve
+ * @return the maximum ordinate value
+ */
+ static uint32_t maxOrdinate(uint32_t level);
+
+ /**
+ * The level of the finite Hilbert curve which contains at least
+ * the given number of points.
+ *
+ * @param numPoints the number of points required
+ * @return the level of the curve
+ */
+ static uint32_t level(uint32_t numPoints);
+
+
+private:
+
+ static uint32_t deinterleave(uint32_t x);
+
+ static uint32_t interleave(uint32_t x);
+
+ static uint32_t prefixScan(uint32_t x);
+
+ static uint32_t descan(uint32_t x);
+
+ static void checkLevel(int level);
+
+
+};
+
+
+} // namespace geos.shape.fractal
+} // namespace geos.shape
+} // namespace geos
+
+
+
diff --git a/include/geos/shape/fractal/HilbertEncoder.h b/include/geos/shape/fractal/HilbertEncoder.h
new file mode 100644
index 0000000..61c0010
--- /dev/null
+++ b/include/geos/shape/fractal/HilbertEncoder.h
@@ -0,0 +1,60 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2020 Paul Ramsey <pramsey at cleverelephant.ca>
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU Lesser General Public Licence as published
+ * by the Free Software Foundation.
+ * See the COPYING file for more information.
+ *
+ **********************************************************************/
+
+
+#pragma once
+
+#include <geos/export.h>
+#include <string>
+#include <vector>
+
+// Forward declarations
+namespace geos {
+namespace geom {
+class Coordinate;
+class Geometry;
+class Envelope;
+}
+}
+
+namespace geos {
+namespace shape { // geos.shape
+namespace fractal { // geos.shape.fractal
+
+
+class GEOS_DLL HilbertEncoder {
+
+public:
+
+ HilbertEncoder(uint32_t p_level, geom::Envelope& extent);
+ uint32_t encode(const geom::Envelope* env);
+ static void sort(std::vector<geom::Geometry*>& geoms);
+
+private:
+
+ uint32_t level;
+ double minx;
+ double miny;
+ double strideX;
+ double strideY;
+
+};
+
+
+} // namespace geos.shape.fractal
+} // namespace geos.shape
+} // namespace geos
+
+
+
diff --git a/include/geos/shape/fractal/Makefile.am b/include/geos/shape/fractal/Makefile.am
new file mode 100644
index 0000000..105a0d3
--- /dev/null
+++ b/include/geos/shape/fractal/Makefile.am
@@ -0,0 +1,13 @@
+#
+# This file is part of project GEOS (http://trac.osgeo.org/geos/)
+#
+SUBDIRS =
+
+EXTRA_DIST =
+
+geosdir = $(includedir)/geos/shape/fractal
+
+geos_HEADERS = \
+ HilbertCode.h \
+ HilbertEncoder.h \
+ MortonCode.h
diff --git a/include/geos/shape/fractal/MortonCode.h b/include/geos/shape/fractal/MortonCode.h
new file mode 100644
index 0000000..6743f87
--- /dev/null
+++ b/include/geos/shape/fractal/MortonCode.h
@@ -0,0 +1,139 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2020 Paul Ramsey <pramsey at cleverelephant.ca>
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU Lesser General Public Licence as published
+ * by the Free Software Foundation.
+ * See the COPYING file for more information.
+ *
+ **********************************************************************/
+
+
+#pragma once
+
+#include <geos/export.h>
+#include <string>
+
+// Forward declarations
+namespace geos {
+namespace geom {
+class Coordinate;
+}
+}
+
+namespace geos {
+namespace shape { // geos.shape
+namespace fractal { // geos.shape.fractal
+
+
+/**
+ * Encodes points as the index along the planar Morton (Z-order) curve.
+ *
+ * The planar Morton (Z-order) curve is a continuous space-filling curve.
+ * The Morton curve defines an ordering of the
+ * points in the positive quadrant of the plane.
+ * The index of a point along the Morton curve is called the Morton code.
+ *
+ * A sequence of subsets of the Morton curve can be defined by a level number.
+ * Each level subset occupies a square range.
+ * The curve at level n M(n) contains 2^(n + 1) points.
+ * It fills the range square of side 2^level.
+ * Curve points have ordinates in the range [0, 2^level - 1].
+ * The code for a given point is identical at all levels.
+ * The level simply determines the number of points in the curve subset
+ * and the size of the range square.
+ *
+ * This implementation represents codes using 32-bit integers.
+ * This allows levels 0 to 16 to be handled.
+ * The class supports encoding points
+ * and decoding the point for a given code value.
+ *
+ * The Morton order has the property that it tends to preserve locality.
+ * This means that codes which are near in value will have spatially proximate
+ * points. The converse is not always true - the delta between
+ * codes for nearby points is not always small. But the average delta
+ * is small enough that the Morton order is an effective way of linearizing space
+ * to support range queries.
+ *
+ * @author Martin Davis
+ *
+ * @see HilbertCode
+ */
+class GEOS_DLL MortonCode {
+
+public:
+
+ /**
+ * The maximum curve level that can be represented.
+ */
+ static constexpr int MAX_LEVEL = 16;
+
+ /**
+ * Computes the index of the point (x,y)
+ * in the Morton curve ordering.
+ *
+ * @param x the x ordinate of the point
+ * @param y the y ordinate of the point
+ * @return the index of the point along the Morton curve
+ */
+ static uint32_t encode(int x, int y);
+
+ /**
+ * Computes the point on the Morton curve
+ * for a given index.
+ *
+ * @param index the index of the point on the curve
+ * @return the point on the curve
+ */
+ static geom::Coordinate decode(uint32_t index);
+
+ /**
+ * The number of points in the curve for the given level.
+ * The number of points is 2<sup>2 * level</sup>.
+ *
+ * @param level the level of the curve
+ * @return the number of points
+ */
+ static uint32_t levelSize(uint32_t level);
+
+ /**
+ * The maximum ordinate value for points
+ * in the curve for the given level.
+ * The maximum ordinate is 2<sup>level</sup> - 1.
+ *
+ * @param level the level of the curve
+ * @return the maximum ordinate value
+ */
+ static uint32_t maxOrdinate(uint32_t level);
+
+ /**
+ * The level of the finite Morton curve which contains at least
+ * the given number of points.
+ *
+ * @param numPoints the number of points required
+ * @return the level of the curve
+ */
+ static uint32_t level(uint32_t numPoints);
+
+
+private:
+
+ static void checkLevel(uint32_t level) ;
+
+ static uint32_t interleave(uint32_t x);
+
+ static uint32_t deinterleave(uint32_t x);
+
+};
+
+
+} // namespace geos.shape.fractal
+} // namespace geos.shape
+} // namespace geos
+
+
+
diff --git a/src/Makefile.am b/src/Makefile.am
index 2e5bc25..38a8b2f 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -16,6 +16,7 @@ SUBDIRS = \
precision \
simplify \
triangulate \
+ shape \
util
EXTRA_DIST = dirlist.mk CMakeLists.txt info.plist.in
@@ -49,4 +50,5 @@ libgeos_la_LIBADD = \
precision/libprecision.la \
simplify/libsimplify.la \
triangulate/libtriangulate.la \
+ shape/libshape.la \
util/libutil.la
diff --git a/src/operation/union/CascadedPolygonUnion.cpp b/src/operation/union/CascadedPolygonUnion.cpp
index 6d1d9bf..5b7361c 100644
--- a/src/operation/union/CascadedPolygonUnion.cpp
+++ b/src/operation/union/CascadedPolygonUnion.cpp
@@ -18,6 +18,7 @@
*
**********************************************************************/
+
#include <geos/operation/union/CascadedPolygonUnion.h>
#include <geos/operation/union/OverlapUnion.h>
#include <geos/operation/overlay/OverlayOp.h>
@@ -29,6 +30,7 @@
#include <geos/geom/MultiPolygon.h>
#include <geos/geom/util/PolygonExtracter.h>
#include <geos/index/strtree/STRtree.h>
+
// std
#include <cassert>
#include <cstddef>
@@ -153,6 +155,7 @@ CascadedPolygonUnion::Union()
* This makes unioning more efficient, since vertices are more likely
* to be eliminated on each round.
*/
+
index::strtree::STRtree index(STRTREE_NODE_CAPACITY);
typedef std::vector<geom::Polygon*>::iterator iterator_type;
diff --git a/src/shape/Makefile.am b/src/shape/Makefile.am
new file mode 100644
index 0000000..8bdb8aa
--- /dev/null
+++ b/src/shape/Makefile.am
@@ -0,0 +1,12 @@
+#
+# This file is part of project GEOS (http://trac.osgeo.org/geos/)
+#
+SUBDIRS = \
+ fractal
+
+noinst_LTLIBRARIES = libshape.la
+
+libshape_la_SOURCES =
+
+libshape_la_LIBADD = \
+ fractal/libshapefractal.la
diff --git a/src/shape/fractal/HilbertCode.cpp b/src/shape/fractal/HilbertCode.cpp
new file mode 100644
index 0000000..484f475
--- /dev/null
+++ b/src/shape/fractal/HilbertCode.cpp
@@ -0,0 +1,199 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2020 Paul Ramsey <pramsey at cleverelephant.ca>
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU Lesser General Public Licence as published
+ * by the Free Software Foundation.
+ * See the COPYING file for more information.
+ *
+ **********************************************************************/
+
+#include <geos/shape/fractal/HilbertCode.h>
+#include <geos/geom/Coordinate.h>
+#include <geos/util/IllegalArgumentException.h>
+
+#include <cstdint>
+
+namespace geos {
+namespace shape { // geos.shape
+namespace fractal { // geos.shape.fractal
+
+// Fast Hilbert curve algorithm by http://threadlocalmutex.com/
+// From C++ https://github.com/rawrunprotected/hilbert_curves
+//(public domain)
+
+/*public static*/
+uint32_t HilbertCode::levelSize(uint32_t level)
+{
+ return (uint32_t)std::pow(2, 2*level);
+}
+
+/*public static*/
+uint32_t HilbertCode::maxOrdinate(uint32_t level)
+{
+ return (uint32_t)std::pow(2, level) - 1;
+}
+
+/*public static*/
+uint32_t HilbertCode::level(uint32_t numPoints)
+{
+ uint32_t pow2 = (uint32_t) ((std::log(numPoints)/std::log(2)));
+ uint32_t level = pow2 / 2;
+ uint32_t size = levelSize(level);
+ if (size < numPoints)
+ level += 1;
+ return level;
+}
+
+
+uint32_t
+HilbertCode::deinterleave(uint32_t x)
+{
+ x = x & 0x55555555;
+ x = (x | (x >> 1)) & 0x33333333;
+ x = (x | (x >> 2)) & 0x0F0F0F0F;
+ x = (x | (x >> 4)) & 0x00FF00FF;
+ x = (x | (x >> 8)) & 0x0000FFFF;
+ return x;
+}
+
+uint32_t
+HilbertCode::interleave(uint32_t x)
+{
+ x = (x | (x << 8)) & 0x00FF00FF;
+ x = (x | (x << 4)) & 0x0F0F0F0F;
+ x = (x | (x << 2)) & 0x33333333;
+ x = (x | (x << 1)) & 0x55555555;
+ return x;
+}
+
+uint32_t
+HilbertCode::prefixScan(uint32_t x)
+{
+ x = (x >> 8) ^ x;
+ x = (x >> 4) ^ x;
+ x = (x >> 2) ^ x;
+ x = (x >> 1) ^ x;
+ return x;
+}
+
+uint32_t
+HilbertCode::descan(uint32_t x)
+{
+ return x ^ (x >> 1);
+}
+
+/*private static*/
+void
+HilbertCode::checkLevel(int level)
+{
+ if (level > MAX_LEVEL) {
+ throw util::IllegalArgumentException("Level out of range");
+ }
+}
+
+/* public static */
+geom::Coordinate
+HilbertCode::decode(uint32_t level, uint32_t i)
+{
+ checkLevel(level);
+ i = i << (32 - 2 * level);
+
+ uint32_t i0 = deinterleave(i);
+ uint32_t i1 = deinterleave(i >> 1);
+
+ uint32_t t0 = (i0 | i1) ^ 0xFFFF;
+ uint32_t t1 = i0 & i1;
+
+ uint32_t prefixT0 = prefixScan(t0);
+ uint32_t prefixT1 = prefixScan(t1);
+
+ uint32_t a = (((i0 ^ 0xFFFF) & prefixT1) | (i0 & prefixT0));
+
+ geom::Coordinate c;
+ c.x = (a ^ i1) >> (16 - level);
+ c.y = (a ^ i0 ^ i1) >> (16 - level);
+ return c;
+}
+
+/* public static */
+uint32_t
+HilbertCode::encode(uint32_t level, uint32_t x, uint32_t y)
+{
+ checkLevel(level);
+ x = x << (16 - level);
+ y = y << (16 - level);
+
+ uint32_t A, B, C, D;
+
+ // Initial prefix scan round, prime with x and y
+ {
+ uint32_t a = x ^ y;
+ uint32_t b = 0xFFFF ^ a;
+ uint32_t c = 0xFFFF ^ (x | y);
+ uint32_t d = x & (y ^ 0xFFFF);
+
+ A = a | (b >> 1);
+ B = (a >> 1) ^ a;
+
+ C = ((c >> 1) ^ (b & (d >> 1))) ^ c;
+ D = ((a & (c >> 1)) ^ (d >> 1)) ^ d;
+ }
+
+ {
+ uint32_t a = A;
+ uint32_t b = B;
+ uint32_t c = C;
+ uint32_t d = D;
+
+ A = ((a & (a >> 2)) ^ (b & (b >> 2)));
+ B = ((a & (b >> 2)) ^ (b & ((a ^ b) >> 2)));
+
+ C ^= ((a & (c >> 2)) ^ (b & (d >> 2)));
+ D ^= ((b & (c >> 2)) ^ ((a ^ b) & (d >> 2)));
+ }
+
+ {
+ uint32_t a = A;
+ uint32_t b = B;
+ uint32_t c = C;
+ uint32_t d = D;
+
+ A = ((a & (a >> 4)) ^ (b & (b >> 4)));
+ B = ((a & (b >> 4)) ^ (b & ((a ^ b) >> 4)));
+
+ C ^= ((a & (c >> 4)) ^ (b & (d >> 4)));
+ D ^= ((b & (c >> 4)) ^ ((a ^ b) & (d >> 4)));
+ }
+
+ // Final round and projection
+ {
+ uint32_t a = A;
+ uint32_t b = B;
+ uint32_t c = C;
+ uint32_t d = D;
+
+ C ^= ((a & (c >> 8)) ^ (b & (d >> 8)));
+ D ^= ((b & (c >> 8)) ^ ((a ^ b) & (d >> 8)));
+ }
+
+ // Undo transformation prefix scan
+ uint32_t a = C ^ (C >> 1);
+ uint32_t b = D ^ (D >> 1);
+
+ // Recover index bits
+ uint32_t i0 = x ^ y;
+ uint32_t i1 = b | (0xFFFF ^ (i0 | a));
+
+ return ((interleave(i1) << 1) | interleave(i0)) >> (32 - 2 * level);
+}
+
+
+
+} // namespace geos.shape.fractal
+} // namespace geos.shape
+} // namespace geos
diff --git a/src/shape/fractal/HilbertEncoder.cpp b/src/shape/fractal/HilbertEncoder.cpp
new file mode 100644
index 0000000..d853b7c
--- /dev/null
+++ b/src/shape/fractal/HilbertEncoder.cpp
@@ -0,0 +1,91 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2020 Paul Ramsey <pramsey at cleverelephant.ca>
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU Lesser General Public Licence as published
+ * by the Free Software Foundation.
+ * See the COPYING file for more information.
+ *
+ **********************************************************************/
+
+#include <geos/shape/fractal/HilbertEncoder.h>
+#include <geos/shape/fractal/HilbertCode.h>
+#include <geos/geom/Coordinate.h>
+#include <geos/geom/Envelope.h>
+#include <geos/geom/Geometry.h>
+
+#include <cstdint>
+#include <vector>
+
+namespace geos {
+namespace shape { // geos.shape
+namespace fractal { // geos.shape.fractal
+
+
+HilbertEncoder::HilbertEncoder(uint32_t p_level, geom::Envelope& extent)
+ : level(p_level)
+{
+ int hside = (int)std::pow(2, level) - 1;
+
+ minx = extent.getMinX();
+ strideX = extent.getWidth() / hside;
+
+ miny = extent.getMinY();
+ strideY = extent.getHeight() / hside;
+}
+
+uint32_t
+HilbertEncoder::encode(const geom::Envelope* env)
+{
+ double midx = env->getWidth()/2 + env->getMinX();
+ uint32_t x = (uint32_t) ((midx - minx) / strideX);
+
+ double midy = env->getHeight()/2 + env->getMinY();
+ uint32_t y = (uint32_t) ((midy - miny) / strideY);
+
+ return HilbertCode::encode(level, x, y);
+}
+
+
+/* public static */
+void
+HilbertEncoder::sort(std::vector<geom::Geometry*>& geoms)
+{
+ struct HilbertComparator {
+
+ HilbertEncoder& enc;
+
+ HilbertComparator(HilbertEncoder& e)
+ : enc(e) {};
+
+ bool
+ operator()(const geom::Geometry* a, const geom::Geometry* b)
+ {
+ return enc.encode(a->getEnvelopeInternal()) > enc.encode(b->getEnvelopeInternal());
+ }
+ };
+
+ geom::Envelope extent;
+ for (const geom::Geometry* geom: geoms)
+ {
+ if (extent.isNull())
+ extent = *(geom->getEnvelopeInternal());
+ else
+ extent.expandToInclude(*(geom->getEnvelopeInternal()));
+ }
+ if (extent.isNull()) return;
+
+ HilbertEncoder encoder(12, extent);
+ HilbertComparator hilbertCompare(encoder);
+ std::sort(geoms.begin(), geoms.end(), hilbertCompare);
+ return;
+}
+
+
+} // namespace geos.shape.fractal
+} // namespace geos.shape
+} // namespace geos
diff --git a/src/shape/fractal/Makefile.am b/src/shape/fractal/Makefile.am
new file mode 100644
index 0000000..4dfa260
--- /dev/null
+++ b/src/shape/fractal/Makefile.am
@@ -0,0 +1,13 @@
+#
+# This file is part of project GEOS (http://trac.osgeo.org/geos/)
+#
+noinst_LTLIBRARIES = libshapefractal.la
+
+AM_CPPFLAGS = -I$(top_srcdir)/include
+
+libshapefractal_la_SOURCES = \
+ HilbertCode.cpp \
+ HilbertEncoder.cpp \
+ MortonCode.cpp
+
+libshapefractal_la_LIBADD =
diff --git a/src/shape/fractal/MortonCode.cpp b/src/shape/fractal/MortonCode.cpp
new file mode 100644
index 0000000..bf1715a
--- /dev/null
+++ b/src/shape/fractal/MortonCode.cpp
@@ -0,0 +1,136 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2020 Paul Ramsey <pramsey at cleverelephant.ca>
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU Lesser General Public Licence as published
+ * by the Free Software Foundation.
+ * See the COPYING file for more information.
+ *
+ **********************************************************************/
+
+#include <geos/shape/fractal/MortonCode.h>
+#include <geos/geom/Coordinate.h>
+#include <geos/util/IllegalArgumentException.h>
+
+#include <cstdint>
+
+namespace geos {
+namespace shape { // geos.shape
+namespace fractal { // geos.shape.fractal
+
+
+/**
+* The number of points in the curve for the given level.
+* The number of points is 2<sup>2 * level</sup>.
+*
+* @param level the level of the curve
+* @return the number of points
+*/
+uint32_t
+MortonCode::levelSize(uint32_t level)
+{
+ checkLevel(level);
+ return (uint32_t) std::pow(2, 2 * level);
+}
+
+/**
+* The maximum ordinate value for points
+* in the curve for the given level.
+* The maximum ordinate is 2<sup>level</sup> - 1.
+*
+* @param level the level of the curve
+* @return the maximum ordinate value
+*/
+uint32_t
+MortonCode::maxOrdinate(uint32_t level)
+{
+ checkLevel(level);
+ return (uint32_t) std::pow(2, level) - 1;
+}
+
+/**
+* The level of the finite Morton curve which contains at least
+* the given number of points.
+*
+* @param numPoints the number of points required
+* @return the level of the curve
+*/
+uint32_t
+MortonCode::level(uint32_t numPoints)
+{
+ uint32_t pow2 = (uint32_t) ( (std::log(numPoints)/std::log(2)));
+ uint32_t level = pow2 / 2;
+ uint32_t sz = levelSize(level);
+ if (sz < numPoints)
+ level += 1;
+ return level;
+}
+
+void
+MortonCode::checkLevel(uint32_t level)
+{
+ if (level > MAX_LEVEL) {
+ throw util::IllegalArgumentException("Level not in range");
+ }
+}
+
+/**
+* Computes the index of the point (x,y)
+* in the Morton curve ordering.
+*
+* @param x the x ordinate of the point
+* @param y the y ordinate of the point
+* @return the index of the point along the Morton curve
+*/
+uint32_t
+MortonCode::encode(int x, int y)
+{
+ return (interleave(y) << 1) + interleave(x);
+}
+
+uint32_t
+MortonCode::interleave(uint32_t x)
+{
+ x &= 0x0000ffff; // x = ---- ---- ---- ---- fedc ba98 7654 3210
+ x = (x ^ (x << 8)) & 0x00ff00ff; // x = ---- ---- fedc ba98 ---- ---- 7654 3210
+ x = (x ^ (x << 4)) & 0x0f0f0f0f; // x = ---- fedc ---- ba98 ---- 7654 ---- 3210
+ x = (x ^ (x << 2)) & 0x33333333; // x = --fe --dc --ba --98 --76 --54 --32 --10
+ x = (x ^ (x << 1)) & 0x55555555; // x = -f-e -d-c -b-a -9-8 -7-6 -5-4 -3-2 -1-0
+ return x;
+}
+
+/**
+* Computes the point on the Morton curve
+* for a given index.
+*
+* @param index the index of the point on the curve
+* @return the point on the curve
+*/
+geom::Coordinate
+MortonCode::decode(uint32_t index)
+{
+ uint32_t x = deinterleave(index);
+ uint32_t y = deinterleave(index >> 1);
+ return geom::Coordinate(x, y);
+}
+
+
+uint32_t
+MortonCode::deinterleave(uint32_t x)
+{
+ x = x & 0x55555555;
+ x = (x | (x >> 1)) & 0x33333333;
+ x = (x | (x >> 2)) & 0x0F0F0F0F;
+ x = (x | (x >> 4)) & 0x00FF00FF;
+ x = (x | (x >> 8)) & 0x0000FFFF;
+ return x;
+}
+
+
+} // namespace geos.shape.fractal
+} // namespace geos.shape
+} // namespace geos
diff --git a/tests/unit/Makefile.am b/tests/unit/Makefile.am
index fa5a2a7..51129c8 100644
--- a/tests/unit/Makefile.am
+++ b/tests/unit/Makefile.am
@@ -226,6 +226,8 @@ geos_unit_SOURCES = \
triangulate/quadedge/QuadEdgeTest.cpp \
triangulate/quadedge/VertexTest.cpp \
triangulate/VoronoiTest.cpp \
+ shape/fractal/HilbertCodeTest.cpp \
+ shape/fractal/MortonCodeTest.cpp \
util/NodingTestUtil.cpp \
util/UniqueCoordinateArrayFilterTest.cpp
diff --git a/tests/unit/shape/fractal/HilbertCodeTest.cpp b/tests/unit/shape/fractal/HilbertCodeTest.cpp
new file mode 100644
index 0000000..b28ebaa
--- /dev/null
+++ b/tests/unit/shape/fractal/HilbertCodeTest.cpp
@@ -0,0 +1,125 @@
+// Test Suite for HilbertCode class.
+
+// tut
+#include <tut/tut.hpp>
+// geos
+#include <geos/constants.h>
+#include <geos/shape/fractal/HilbertCode.h>
+#include <geos/geom/Coordinate.h>
+// std
+
+using geos::shape::fractal::HilbertCode;
+using geos::geom::Coordinate;
+
+namespace tut {
+
+// Common data used by tests
+struct test_hilbertcode_data {
+
+ void checkDecode(uint32_t order, uint32_t index, int x, int y)
+ {
+ Coordinate p = HilbertCode::decode(order, index);
+ //System.out.println(p);
+ ensure_equals(x, (int) p.x);
+ ensure_equals(y, (int) p.y);
+ }
+
+ void checkDecodeEncodeForLevel(uint32_t level)
+ {
+ uint32_t n = HilbertCode::levelSize(level);
+ for (int i = 0; i < n; i++) {
+ checkDecodeEncode(level, i);
+ }
+ }
+
+
+ void
+ checkDecodeEncode(uint32_t level, uint32_t index)
+ {
+ Coordinate p = HilbertCode::decode(level, index);
+ uint32_t encode = HilbertCode::encode(level, (uint32_t)(p.x), (uint32_t)(p.y));
+ ensure_equals(index, encode);
+ }
+
+};
+
+typedef test_group<test_hilbertcode_data> group;
+typedef group::object object;
+
+group test_hilbertcode_group("geos::shape::fractal::HilbertCode");
+
+template<>
+template<>
+void object::test<1>
+()
+{
+ ensure_equals( HilbertCode::levelSize( 0 ), 1);
+ ensure_equals( HilbertCode::levelSize( 1 ), 4);
+ ensure_equals( HilbertCode::levelSize( 2 ), 16);
+ ensure_equals( HilbertCode::levelSize( 3 ), 64);
+ ensure_equals( HilbertCode::levelSize( 4 ), 256);
+ ensure_equals( HilbertCode::levelSize( 5 ), 1024);
+ ensure_equals( HilbertCode::levelSize( 6 ), 4096);
+}
+
+
+template<>
+template<>
+void object::test<2>
+()
+{
+ ensure_equals( HilbertCode::level( 1 ), 0);
+
+ ensure_equals( HilbertCode::level( 2 ), 1);
+ ensure_equals( HilbertCode::level( 3 ), 1);
+ ensure_equals( HilbertCode::level( 4 ), 1);
+
+ ensure_equals( HilbertCode::level( 5 ), 2);
+ ensure_equals( HilbertCode::level( 13 ), 2);
+ ensure_equals( HilbertCode::level( 15 ), 2);
+ ensure_equals( HilbertCode::level( 16 ), 2);
+
+ ensure_equals( HilbertCode::level( 17 ), 3);
+ ensure_equals( HilbertCode::level( 63 ), 3);
+ ensure_equals( HilbertCode::level( 64 ), 3);
+
+ ensure_equals( HilbertCode::level( 65 ), 4);
+ ensure_equals( HilbertCode::level( 255 ), 4);
+ ensure_equals( HilbertCode::level( 255 ), 4);
+ ensure_equals( HilbertCode::level( 256 ), 4);
+}
+
+
+template<>
+template<>
+void object::test<3>
+()
+{
+ checkDecode(1, 0, 0, 0);
+
+ checkDecode(1, 0, 0, 0);
+ checkDecode(1, 1, 0, 1);
+
+ checkDecode(3, 0, 0, 0);
+ checkDecode(3, 1, 0, 1);
+
+ checkDecode(4,0, 0, 0);
+ checkDecode(4, 1, 1, 0);
+ checkDecode(4, 24, 6, 2);
+ checkDecode(4, 255, 15, 0);
+
+ checkDecode(5, 124, 8, 6);
+}
+
+
+template<>
+template<>
+void object::test<4>
+()
+{
+ checkDecodeEncodeForLevel(4);
+ checkDecodeEncodeForLevel(5);
+}
+
+
+} // namespace tut
diff --git a/tests/unit/shape/fractal/MortonCodeTest.cpp b/tests/unit/shape/fractal/MortonCodeTest.cpp
new file mode 100644
index 0000000..b606653
--- /dev/null
+++ b/tests/unit/shape/fractal/MortonCodeTest.cpp
@@ -0,0 +1,119 @@
+// Test Suite for MortonCode class.
+
+// tut
+#include <tut/tut.hpp>
+// geos
+#include <geos/constants.h>
+#include <geos/shape/fractal/MortonCode.h>
+#include <geos/geom/Coordinate.h>
+// std
+
+using geos::shape::fractal::MortonCode;
+using geos::geom::Coordinate;
+
+namespace tut {
+
+// Common data used by tests
+struct test_mortoncode_data {
+
+ void checkDecode(uint32_t index, int x, int y)
+ {
+ Coordinate p = MortonCode::decode(index);
+ //System.out.println(p);
+ ensure_equals(x, (int) p.x);
+ ensure_equals(y, (int) p.y);
+ }
+
+ void checkDecodeEncodeForLevel(uint32_t level)
+ {
+ uint32_t n = MortonCode::levelSize(level);
+ for (int i = 0; i < n; i++) {
+ checkDecodeEncode(i);
+ }
+ }
+
+
+ void
+ checkDecodeEncode(uint32_t index)
+ {
+ Coordinate p = MortonCode::decode(index);
+ uint32_t encode = MortonCode::encode((uint32_t)(p.x), (uint32_t)(p.y));
+ ensure_equals(index, encode);
+ }
+
+};
+
+typedef test_group<test_mortoncode_data> group;
+typedef group::object object;
+
+group test_mortoncode_group("geos::shape::fractal::MortonCode");
+
+template<>
+template<>
+void object::test<1>
+()
+{
+ ensure_equals( MortonCode::levelSize( 0 ), 1);
+ ensure_equals( MortonCode::levelSize( 1 ), 4);
+ ensure_equals( MortonCode::levelSize( 2 ), 16);
+ ensure_equals( MortonCode::levelSize( 3 ), 64);
+ ensure_equals( MortonCode::levelSize( 4 ), 256);
+ ensure_equals( MortonCode::levelSize( 5 ), 1024);
+ ensure_equals( MortonCode::levelSize( 6 ), 4096);
+}
+
+
+template<>
+template<>
+void object::test<2>
+()
+{
+ ensure_equals( MortonCode::level( 1 ), 0);
+
+ ensure_equals( MortonCode::level( 2 ), 1);
+ ensure_equals( MortonCode::level( 3 ), 1);
+ ensure_equals( MortonCode::level( 4 ), 1);
+
+ ensure_equals( MortonCode::level( 5 ), 2);
+ ensure_equals( MortonCode::level( 13 ), 2);
+ ensure_equals( MortonCode::level( 15 ), 2);
+ ensure_equals( MortonCode::level( 16 ), 2);
+
+ ensure_equals( MortonCode::level( 17 ), 3);
+ ensure_equals( MortonCode::level( 63 ), 3);
+ ensure_equals( MortonCode::level( 64 ), 3);
+
+ ensure_equals( MortonCode::level( 65 ), 4);
+ ensure_equals( MortonCode::level( 255 ), 4);
+ ensure_equals( MortonCode::level( 255 ), 4);
+ ensure_equals( MortonCode::level( 256 ), 4);
+}
+
+
+template<>
+template<>
+void object::test<3>
+()
+{
+ checkDecode(0, 0, 0);
+ checkDecode(1, 1, 0);
+ checkDecode(2, 0, 1);
+ checkDecode(3, 1, 1);
+ checkDecode(4, 2, 0);
+
+ checkDecode(24, 4, 2);
+ checkDecode(124, 14, 6);
+ checkDecode(255, 15, 15);
+}
+
+template<>
+template<>
+void object::test<4>
+()
+{
+ checkDecodeEncodeForLevel(4);
+ checkDecodeEncodeForLevel(5);
+}
+
+
+} // namespace tut
-----------------------------------------------------------------------
Summary of changes:
configure.ac | 15 +-
include/geos/Makefile.am | 1 +
include/geos/{math => shape}/Makefile.am | 8 +-
include/geos/shape/fractal/HilbertCode.h | 127 +++++++++++++
.../fractal/HilbertEncoder.h} | 39 ++--
include/geos/shape/fractal/Makefile.am | 13 ++
include/geos/shape/fractal/MortonCode.h | 139 ++++++++++++++
src/Makefile.am | 2 +
src/operation/union/CascadedPolygonUnion.cpp | 3 +
src/shape/Makefile.am | 12 ++
src/shape/fractal/HilbertCode.cpp | 199 +++++++++++++++++++++
src/shape/fractal/HilbertEncoder.cpp | 91 ++++++++++
src/shape/fractal/Makefile.am | 13 ++
src/shape/fractal/MortonCode.cpp | 136 ++++++++++++++
tests/unit/Makefile.am | 2 +
tests/unit/shape/fractal/HilbertCodeTest.cpp | 125 +++++++++++++
tests/unit/shape/fractal/MortonCodeTest.cpp | 119 ++++++++++++
17 files changed, 1028 insertions(+), 16 deletions(-)
copy include/geos/{math => shape}/Makefile.am (53%)
create mode 100644 include/geos/shape/fractal/HilbertCode.h
copy include/geos/{algorithm/PointInRing.h => shape/fractal/HilbertEncoder.h} (50%)
create mode 100644 include/geos/shape/fractal/Makefile.am
create mode 100644 include/geos/shape/fractal/MortonCode.h
create mode 100644 src/shape/Makefile.am
create mode 100644 src/shape/fractal/HilbertCode.cpp
create mode 100644 src/shape/fractal/HilbertEncoder.cpp
create mode 100644 src/shape/fractal/Makefile.am
create mode 100644 src/shape/fractal/MortonCode.cpp
create mode 100644 tests/unit/shape/fractal/HilbertCodeTest.cpp
create mode 100644 tests/unit/shape/fractal/MortonCodeTest.cpp
hooks/post-receive
--
GEOS
More information about the geos-commits
mailing list