Hilbertize an image

Question

For a computer vision app I want to do a mapping of an image, in such a way that every pixel fit hilbert curve, instead of conventional layout. So task could be as follows:

---

Task description

Given square 2D image A with side 2^^N, (N>0), with common Z-order, or row-major representation of pixels in RAM. I need to run a hilbert curve through this image (Like a digital camera, where photodiodes are soldered in a hilbert curve and numbered and accessed in a 1-d ordered manner). Then pixels along this curve just layed down as a regular scanline in image B. So now any consequent samples in output array B now should represent spatially close points in input array A.

---

Example for N=2

Input: array A (imagination of 2D image)
0 1 E F  
3 2 D C     
4 7 8 B    
5 6 9 A 

Input: array A (1D actual layout in memory)
0 1 E F 3 2 D C 4 7 8 B 5 6 9 A

Output: array B (answer)
0 1 2 3 4 5 6 7 8 9 A B C D E F

Note, that values inside arrays are "values". In this example, for array B accidentally they became same as indices. So basically I want index mapping function B[i]=A[?], e.g B2 = A[E]. Hope code golfing helps.

---

Example for N=9:

Taken popular 512*512 test image. Each pixel was threated as 32bit integer, in RGBA format and algorithm proceeded:

Input:

Output:

---

Reverse question:

Here is reversed-task question on codegolf.

In terms of that question, I want to do opposite transform: "I want to Reravel image A, then Unravel to image B". In terms of examples of that question: I want to take photo the lion on the camera with hilbert-curved-sensor (look Output array for lion photo) and get that long "scanline of spatially close pixels" (Input for lion photo, in that question)

---

Reference implementation:

If i take d2xy() function from LGPL sources found here I may write a following program, which works, but maybe not its best way:

#include "hilbert_curve.h"
int main(){
  int *A,*B,i,x,y;
  for(i=0;i<(1<<(N<<1));++i){
    d2xy(N,i,&x,&y);
    B[i]=A[y*(1<<N)+x];
  }
}

---

Resources:

• Spatial indexing with Quadtrees and Hilbert Curves
• hilbert function, which returns a rendering of Hilbert Curves

Answer 1

Jelly, 38 bytes

U1¦U;$;;UN$
Ll2’ð⁽Ƭb3’s2Ç⁸¡⁵DW¤;+\Ḋœị

A monadic Link accepting a square matrix of side \$2^n\$, \$n\$ a positive integer which yields an array of the inputs elements in Hilbert curve order.

Try it online!

How?

First we create [row, column] index deltas describing the U shape with an initial right move (i.e into the \$n=1\$ matrix's top-left corner from the left) - that is [[0, 1], [1, 0], [0, 1], [-1, 0]].

Then we repeatedly perform a transform keeping the current list of deltas to describe the bottom-right quarter of the next and adding the other three quarters as copies with reversed deltas (swapping row changes with column ones where necessary and negating directions of the new top-right quarter).

Finally we prepend a starting index of [1,0] (to the left of the top-left of our [1-indexed] matrix), cumulatively reduce by addition to get the indexes, discard the [1, 0] and index into the input.

How?

U1¦U;$;;UN$ - Link 1, transform deltas to next size: list of current deltas
U1¦         - reverse the first delta (makes bottom-left)
   U;$      - reverse each of those and concatenate with that (adds top-left)
      ;     - concatenate with the input (adds bottom-right)
        UN$ - reverse each of the input deltas and negate (makes top-right)
       ;    - concatenate (adds top-right)

Ll2’ð⁽Ƭb3’s2Ç⁸¡⁵DW¤;+\Ḋœị - Main Link: square matrix (of side 2^n), M
L                          - length of M
 l2                        - log base 2 (=n)
   ’                       - decrement
    ð                      - start a new dyadic chain (i.e. f(n-1, M)) 
     ⁽Ƭ                   - compressed literal = 4258
        b3                 - to base 3 = [1,2,2,1,1,2,0,1]
          ’                - decrement = [0,1,1,0,0,1,-1,0]
           s2              - split into twos = [[0,1],[1,0],[0,1],[-1,0]]
               ¡           - apply repeatedly...
              ⁸            - ...number of times: chain's left argument, n-1
             Ç             - ...what?: last Link as a monad
                           -           - gets all deltas
                   ¤       - nilad followed by link(s) as a nilad:
                ⁵          -   literal ten
                 D         -   to decimal digits = [1,0]
                  W        -   wrap = [[1,0]]
                    ;      - concatenate with the deltas
                     +\    - cumulative reduce using addition (N.B. Ä wont work here) 
                       Ḋ   - dequeue (drop the [1, 0] leading index)
                        œị - multi-dimensional index into M