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