Native support for MOC operations in the Java HEALPix package

[Martin Reinecke]

Hi François-Xavier!

> Very nice result :)
> 
> If I am correct the NESTED ordering uses a Z-order (or Morton) curve (which rely 
> on bit interleaving).

Correct.

> It explains why computing a HEALPix index is very fast and does not depends on 
> the NORDER.

Well - to a certain extent it does depend on the order.
It is true that a few modern CPUs could do the bit interleaving with a
single assembler operation, but the current approach is to use look-up
tables, which deal with 8 orders in one go. So there is still a loop,
but it has maxorder/8 iterations only.

> Do you have the same kind of algorithm for the Peano-Hilbert?

No. The look-up tables would be too large in that case.

> If yes, would you provide a pointer on it?
> The only ones I am aware of use recursion or a 'for' loop from 0 to NORDER...

I'm using these as well, but with some optimizations the conversions
don't take too long.

> If not, what about performances when computing 100_000_000 Peano-Hilbert HEALPix 
> indexes?

On a single core of an Intel Core i3 2120 (several years old by now) I
get around 30 million conversions per second between NEST and Peano
indices for a maximum order of 13, and around 20 million per second for
a maximum order of 29 (measured with the C++ library, I don't have
measurements for Java yet). This is comparable to the cost of an
ang2pix() and probably not a limiting factor in most circumstances.
In the case of MOCs this cost is mostly irrelevant, since it's only
incurred when the MOC is first created from data. All set operations on
MOCs don't need any conversions.

> If performances using Morton and Peano are similar it is (in addition to MOCs) 
> interesting to index data (in databases, ...)

Performance is not dramatically lower, at least.

> since Peano has "better locality-preserving behaviour" than Morton (quoting 
> Wikipedia but I heard/read this elsewhere),

A quantitative paper on this is for example "Analysis of the clustering
properties of the Hilbert space-filling curve" by Moon et al.
Unfortunately I don't have the link at hand.

> so needs less disk access, ...
> 
> Even if it will not become a standard, I think it is worth having it in the Java 
> HEALPix libraries :)

If you want to have a look at the current implementation, I have
appended it below; it needs to be pasted somewhere into the HealpixUtils
class and is surprisingly small :)

Cheers,
  Martin


  static private long nest_peano_helper (long pix, int order, int dir)
    {
    final byte peano_subpix[][][] =
      { { {0,1,3,2}, {3,0,2,1}, {2,3,1,0}, {1,2,0,3},
          {0,3,1,2}, {1,0,2,3}, {2,1,3,0}, {3,2,0,1} },
        { {0,1,3,2}, {1,3,2,0}, {3,2,0,1}, {2,0,1,3},
          {0,2,3,1}, {1,0,2,3}, {3,1,0,2}, {2,3,1,0} } };
    final byte peano_subpath[][][] =
      { { {4,0,6,0}, {7,5,1,1}, {2,4,2,6}, {3,3,7,5},
          {0,2,4,4}, {5,1,5,3}, {6,6,0,2}, {1,7,3,7} },
        { {4,0,0,6}, {5,1,1,7}, {6,2,2,4}, {7,3,3,5},
          {0,4,4,2}, {1,5,5,3}, {2,6,6,0}, {3,7,7,1} } };
    final byte peano_face2path[][] =
      { { 2,5,2,5,3,6,3,6,2,3,2,3 }, { 2,6,2,3,3,5,2,6,2,3,3,5 } };
    final byte peano_face2face[][] =
      { { 0,5,6,11,10,1,4,7,2,3,8,9 }, { 0,5,8,9,6,1,2,7,10,11,4,3 } };

    int face = (int)(pix>>>(2*order));
    long result = 0L;
    byte path = peano_face2path[dir][face];

    for (int shift=2*order-2; shift>=0; shift-=2)
      {
      byte spix = (byte)((pix>>shift) & 0x3);
      result <<= 2;
      result |= peano_subpix[dir][path][spix];
      path=peano_subpath[dir][path][spix];
      }

    return result + (((long)peano_face2face[dir][face])<<(2*order));
    }
  static public long nest2peano(long pix, int order)
    { return nest_peano_helper(pix,order,0); }
  static public long peano2nest(long pix, int order)
    { return nest_peano_helper(pix,order,1); }