Line data    Source code

       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         440 : static uint32_t GDALHilbertCode(uint32_t x, uint32_t y)
      20             : {
      21             :     // Based on public domain code at
      22             :     // https://github.com/rawrunprotected/hilbert_curves
      23             : 
      24         440 :     uint32_t a = x ^ y;
      25         440 :     uint32_t b = 0xFFFF ^ a;
      26         440 :     uint32_t c = 0xFFFF ^ (x | y);
      27         440 :     uint32_t d = x & (y ^ 0xFFFF);
      28             : 
      29         440 :     uint32_t A = a | (b >> 1);
      30         440 :     uint32_t B = (a >> 1) ^ a;
      31         440 :     uint32_t C = ((c >> 1) ^ (b & (d >> 1))) ^ c;
      32         440 :     uint32_t D = ((a & (c >> 1)) ^ (d >> 1)) ^ d;
      33             : 
      34         440 :     a = A;
      35         440 :     b = B;
      36         440 :     c = C;
      37         440 :     d = D;
      38         440 :     A = ((a & (a >> 2)) ^ (b & (b >> 2)));
      39         440 :     B = ((a & (b >> 2)) ^ (b & ((a ^ b) >> 2)));
      40         440 :     C ^= ((a & (c >> 2)) ^ (b & (d >> 2)));
      41         440 :     D ^= ((b & (c >> 2)) ^ ((a ^ b) & (d >> 2)));
      42             : 
      43         440 :     a = A;
      44         440 :     b = B;
      45         440 :     c = C;
      46         440 :     d = D;
      47         440 :     A = ((a & (a >> 4)) ^ (b & (b >> 4)));
      48         440 :     B = ((a & (b >> 4)) ^ (b & ((a ^ b) >> 4)));
      49         440 :     C ^= ((a & (c >> 4)) ^ (b & (d >> 4)));
      50         440 :     D ^= ((b & (c >> 4)) ^ ((a ^ b) & (d >> 4)));
      51             : 
      52         440 :     a = A;
      53         440 :     b = B;
      54         440 :     c = C;
      55         440 :     d = D;
      56         440 :     C ^= ((a & (c >> 8)) ^ (b & (d >> 8)));
      57         440 :     D ^= ((b & (c >> 8)) ^ ((a ^ b) & (d >> 8)));
      58             : 
      59         440 :     a = C ^ (C >> 1);
      60         440 :     b = D ^ (D >> 1);
      61             : 
      62         440 :     uint32_t i0 = x ^ y;
      63         440 :     uint32_t i1 = b | (0xFFFF ^ (i0 | a));
      64             : 
      65         440 :     i0 = (i0 | (i0 << 8)) & 0x00FF00FF;
      66         440 :     i0 = (i0 | (i0 << 4)) & 0x0F0F0F0F;
      67         440 :     i0 = (i0 | (i0 << 2)) & 0x33333333;
      68         440 :     i0 = (i0 | (i0 << 1)) & 0x55555555;
      69             : 
      70         440 :     i1 = (i1 | (i1 << 8)) & 0x00FF00FF;
      71         440 :     i1 = (i1 | (i1 << 4)) & 0x0F0F0F0F;
      72         440 :     i1 = (i1 | (i1 << 2)) & 0x33333333;
      73         440 :     i1 = (i1 | (i1 << 1)) & 0x55555555;
      74             : 
      75         440 :     uint32_t value = ((i1 << 1) | i0);
      76             : 
      77         440 :     return value;
      78             : }
      79             : 
      80         440 : uint32_t GDALHilbertCode(const OGREnvelope *poDomain, double dfX, double dfY)
      81             : {
      82         440 :     uint32_t x = 0;
      83         440 :     uint32_t y = 0;
      84         440 :     if (poDomain->Width() != 0.0)
      85         440 :         x = static_cast<uint32_t>(std::round(
      86         440 :             GDAL_HILBERT_MAX * (dfX - poDomain->MinX) / poDomain->Width()));
      87         440 :     if (poDomain->Height() != 0.0)
      88         440 :         y = static_cast<uint32_t>(std::round(
      89         440 :             GDAL_HILBERT_MAX * (dfY - poDomain->MinY) / poDomain->Height()));
      90         440 :     return GDALHilbertCode(x, y);
      91             : }