/**********************************************************************

 *

 * GEOS - Geometry Engine Open Source

 * http://geos.osgeo.org

 *

 * Copyright (C) 2020 Paul Ramsey <pramsey@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(uint32_t 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