/src/gdal/alg/hilbert.cpp
Line | Count | Source |
1 | | /****************************************************************************** |
2 | | * |
3 | | * Project: GDAL Algorithms |
4 | | * Purpose: Hilbert encoding |
5 | | * Author: Björn Harrtell <bjorn at wololo dot org> |
6 | | * |
7 | | ****************************************************************************** |
8 | | * Copyright (c) 2018-2020, Björn Harrtell <bjorn at wololo dot org> |
9 | | * |
10 | | * SPDX-License-Identifier: MIT |
11 | | ****************************************************************************/ |
12 | | |
13 | | #include "gdal_alg.h" |
14 | | |
15 | | // Subtract 1 from the theoretical max, reserving that number for empty or |
16 | | // null geometries. |
17 | | constexpr uint32_t GDAL_HILBERT_MAX = (1 << 16) - 2; |
18 | | |
19 | | static uint32_t GDALHilbertCode(uint32_t x, uint32_t y) |
20 | 0 | { |
21 | | // Based on public domain code at |
22 | | // https://github.com/rawrunprotected/hilbert_curves |
23 | |
|
24 | 0 | uint32_t a = x ^ y; |
25 | 0 | uint32_t b = 0xFFFF ^ a; |
26 | 0 | uint32_t c = 0xFFFF ^ (x | y); |
27 | 0 | uint32_t d = x & (y ^ 0xFFFF); |
28 | |
|
29 | 0 | uint32_t A = a | (b >> 1); |
30 | 0 | uint32_t B = (a >> 1) ^ a; |
31 | 0 | uint32_t C = ((c >> 1) ^ (b & (d >> 1))) ^ c; |
32 | 0 | uint32_t D = ((a & (c >> 1)) ^ (d >> 1)) ^ d; |
33 | |
|
34 | 0 | a = A; |
35 | 0 | b = B; |
36 | 0 | c = C; |
37 | 0 | d = D; |
38 | 0 | A = ((a & (a >> 2)) ^ (b & (b >> 2))); |
39 | 0 | B = ((a & (b >> 2)) ^ (b & ((a ^ b) >> 2))); |
40 | 0 | C ^= ((a & (c >> 2)) ^ (b & (d >> 2))); |
41 | 0 | D ^= ((b & (c >> 2)) ^ ((a ^ b) & (d >> 2))); |
42 | |
|
43 | 0 | a = A; |
44 | 0 | b = B; |
45 | 0 | c = C; |
46 | 0 | d = D; |
47 | 0 | A = ((a & (a >> 4)) ^ (b & (b >> 4))); |
48 | 0 | B = ((a & (b >> 4)) ^ (b & ((a ^ b) >> 4))); |
49 | 0 | C ^= ((a & (c >> 4)) ^ (b & (d >> 4))); |
50 | 0 | D ^= ((b & (c >> 4)) ^ ((a ^ b) & (d >> 4))); |
51 | |
|
52 | 0 | a = A; |
53 | 0 | b = B; |
54 | 0 | c = C; |
55 | 0 | d = D; |
56 | 0 | C ^= ((a & (c >> 8)) ^ (b & (d >> 8))); |
57 | 0 | D ^= ((b & (c >> 8)) ^ ((a ^ b) & (d >> 8))); |
58 | |
|
59 | 0 | a = C ^ (C >> 1); |
60 | 0 | b = D ^ (D >> 1); |
61 | |
|
62 | 0 | uint32_t i0 = x ^ y; |
63 | 0 | uint32_t i1 = b | (0xFFFF ^ (i0 | a)); |
64 | |
|
65 | 0 | i0 = (i0 | (i0 << 8)) & 0x00FF00FF; |
66 | 0 | i0 = (i0 | (i0 << 4)) & 0x0F0F0F0F; |
67 | 0 | i0 = (i0 | (i0 << 2)) & 0x33333333; |
68 | 0 | i0 = (i0 | (i0 << 1)) & 0x55555555; |
69 | |
|
70 | 0 | i1 = (i1 | (i1 << 8)) & 0x00FF00FF; |
71 | 0 | i1 = (i1 | (i1 << 4)) & 0x0F0F0F0F; |
72 | 0 | i1 = (i1 | (i1 << 2)) & 0x33333333; |
73 | 0 | i1 = (i1 | (i1 << 1)) & 0x55555555; |
74 | |
|
75 | 0 | uint32_t value = ((i1 << 1) | i0); |
76 | |
|
77 | 0 | return value; |
78 | 0 | } |
79 | | |
80 | | uint32_t GDALHilbertCode(const OGREnvelope *poDomain, double dfX, double dfY) |
81 | 0 | { |
82 | 0 | uint32_t x = 0; |
83 | 0 | uint32_t y = 0; |
84 | 0 | if (poDomain->Width() != 0.0) |
85 | 0 | x = static_cast<uint32_t>(std::round( |
86 | 0 | GDAL_HILBERT_MAX * (dfX - poDomain->MinX) / poDomain->Width())); |
87 | 0 | if (poDomain->Height() != 0.0) |
88 | 0 | y = static_cast<uint32_t>(std::round( |
89 | 0 | GDAL_HILBERT_MAX * (dfY - poDomain->MinY) / poDomain->Height())); |
90 | 0 | return GDALHilbertCode(x, y); |
91 | 0 | } |