/* 
   ======================================================================
   Routines to generate and compare Morton Codes, i.e. indices on space
   filling Z-curves.
   
   Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) 
   (c) Lucerne University of Applied Sciences and Arts,
       supported by the Swiss National Science Foundation 
       (SNSF #147053, "Daylight Redirecting Components",
        SNSF #179067, "Light Fields for Spatio-Temporal Glare Assessment")
   ======================================================================   
   
   $Id$
*/

#ifndef _MORTON_H
   #define _MORTON_H
   
   #include "fvect.h"
   #include <stdint.h>

   
   /* Type to represent Morton codes (Z-curve indices) */
   typedef uint64_t MortonIdx;
   
   
   /* Number of bits per dimension for Morton code (Z-curve index) in 2D 
      and 3D. This corresponds to the maximum number of bounding box
      subdivisions which can be resolved per axis; further subdivisions 
      will map to the same Morton code. */
   #define MORTON3D_BITS 21
   #define MORTON3D_MAX  (((MortonIdx)1 << MORTON3D_BITS) - 1)

   #define MORTON2D_BITS 32
   #define MORTON2D_MAX  (((MortonIdx)1 << MORTON2D_BITS) - 1)


   /* Interleave lower MORTONxx_BITS bits of k with 00 for 3D (resp. 0 for 
      2D), resulting in 3*MORTON3D_BITS (resp. 2*MORTON2D_BITS) bits. This
      is an optimised "bitmask hack" version code taken from: 
         http://www.forceflow.be/2013/10/07/
            morton-encodingdecoding-through-bit-interleaving-implementations/
      and (for 32-bit codes):
         https://fgiesen.wordpress.com/2009/12/13/decoding-morton-codes/
   */
   #define Morton3DBitInterleave(k) ((k) &= MORTON3D_MAX, \
      (k) = ((k) | (k) << 32) & 0x001f00000000ffff, \
      (k) = ((k) | (k) << 16) & 0x001f0000ff0000ff, \
      (k) = ((k) | (k) << 8)  & 0x100f00f00f00f00f, \
      (k) = ((k) | (k) << 4)  & 0x10c30c30c30c30c3, \
      (k) = ((k) | (k) << 2)  & 0x1249249249249249  \
   )

   #define Morton2DBitInterleave(k) ((k) &= MORTON2D_MAX, \
      (k) = ((k) | (k) << 16) & 0x0000ffff0000ffff, \
      (k) = ((k) | (k) << 8)  & 0x00ff00ff00ff00ff, \
      (k) = ((k) | (k) << 4)  & 0x0f0f0f0f0f0f0f0f, \
      (k) = ((k) | (k) << 2)  & 0x3333333333333333, \
      (k) = ((k) | (k) << 1)  & 0x5555555555555555  \
   )


   /* Compute Morton code (Z-curve index) of length 3 * MORTON3D_BITS bits 
      for 3D key within bounding box defined by org and scaled to maximum
      index with scale */
   MortonIdx Key2Morton3D (const FVECT key, const FVECT org, RREAL scale);
   
   /* As above, but compute Morton code of length 2 * MORTON2D_BITS bits 
      for 2D key */
   MortonIdx Key2Morton2D (const RREAL key [2], const RREAL org [2], 
      RREAL scale
   );
#endif /* _MORTON_H */