Hilbert-Curvify a Matrix Inspired by this question

Another way to unroll a 2D image into a 1D string is to use an Hilbert Curve.

There are many version of this curve, depending on the number of iterations used while computing it. Below follow example of Hilbert Curves from first order to fifth order.



The way of computing this curve is the following. First we define the first order Hilbert Curve as the one shown in figure (the one for n = 1), so that it fits in a 1x1 square. We than make four copies of this curve, spacing them in a 4x4 square, so that they all present the "concavity" towards the left side. We then flip the two leftmost order 1 curves, so that the top one concavity faces toward the top, while the bottom's faces the bottom. We finally connect the the corners of the adjacent Hilbert Curves. If wanting to obtain a (n+1)-order Curve, we just need to repeat the process with four n-order Curves. We can see a visualisation of the process here (I will also add an image detailing the process soon)

Your task in this challenge is to unroll a matrix of integers along the lowest order Hilbert Curve for that matrix.

For simplicity's sake, we will have the curve starting from the top left corner of the matrix.

You can receive the input either as a list of list of integers, where each
sub-list represents a row of the matrix.

You can assume that the input will be a square matrix (n*n).

For example:

Input:

[[ 1, 2,]
 [ 3, 4 ]]


Output:

[ 1, 2, 4, 3 ]


Since we are using the first order Hilbert Curve shown in figure

Input:

[[ 1, 2, 3, 4,    ]
 [ 5, 6, 7, 8,    ]
 [ 9, 10, 11, 12, ]
 [ 13, 14, 15, 16 ]]


Output:

[ 1, 5, 6, 2, 3, 4, 8, 7, 11, 12, 16, 15, 14, 10, 9, 13 ]


Using the second order Hilbert Curve

As usual, standard loopholes are not permitted.

This is code-golf, so the shortest answer in byte wins.
APL (Dyalog Unicode), 41 bytesSBCS

Saved 30 bytes (!) by consulting the wisdom of the APL Orchard, especially @ngn and @Sherlock9.


{0::⍵⋄∊∇¨⌽∘⊖¨@4,⌽@1⊢∘⍉\⌽↑∘⍵¨∘.,⍨2 ¯2÷⍨≢⍵}


Try it online!

Explanation as follows:

{0::⍵⋄∊∇¨⌽∘⊖¨@4,⌽@1⊢∘⍉\⌽↑∘⍵¨∘.,⍨2 ¯2÷⍨≢⍵} ⍝ Recursive function - takes input as an
                                          ⍝ n*n square matrix
 0::⍵                                     ⍝ Our base case - this is an error guard
                                          ⍝ If there's any error, catch it and
                                          ⍝ return the function's input
                                      ≢⍵  ⍝ Find the number of rows in the input
                                2 ¯2÷⍨    ⍝ Divide the above by 2 and negative 2,
                                          ⍝ resulting in a 2-element vector
                            ∘.,⍨          ⍝ Outer product - take the above vector and
                                          ⍝ apply concatenation (,) with each element
                                          ⍝ against all elements in the vector. Since
                                          ⍝ we have a 2-element vector, this results in
                                          ⍝ a 2-by-2 matrix, e.g.
                                          ⍝ [[(2,2),(2,¯2)],[(¯2,2),(¯2,¯2)]]
                        ↑∘⍵¨              ⍝ For each element in the matrix, we apply
                                          ⍝ "take" against our original input matrix.
                                          ⍝ Take, given a negative number, will take
                                          ⍝ elements from the end of a particular rank.
                                          ⍝ With our argument above, this means that we end
                                          ⍝ up with our original original input matrix
                                          ⍝ split by quadrant into a 2-by-2 matrix.
                                          ⍝ It is also worth noting that take expects
                                          ⍝ an integer argument, so for matrices whose
                                          ⍝ rowcount divided by two results in a decimal
                                          ⍝ (i.e., 1-by-1 matrices), we throw an error
                                          ⍝ which is caught by the guard above, returning
                                          ⍝ the original input.
                       ⌽                  ⍝ Flip the above matrix about the vertical axis.
                   ⊢∘⍉\                   ⍝ Apply a "monadic transpose scan". More details
                                          ⍝ on how this works below, but for our purposes
                                          ⍝ this applies transpose to each of the two 
                                          ⍝ sub-matrices on the right half.
                ⌽@1                       ⍝ Swap the two upper sub-matrices. Given our
                                          ⍝ flip for the overall matrix above, this returns
                                          ⍝ the two upper quadrants to their original
                                          ⍝ positions.
               ,                          ⍝ Ravel: flatten the 2-by-2 matrix into a
                                          ⍝ 4-element vector
         ⌽∘⊖¨@4                           ⍝ Take the last element of the list (the lower
                                          ⍝ right quadrant originally) and flip it
                                          ⍝ along the vertical and horizontal axes. Given
                                          ⍝ the transposition above, this has the final
                                          ⍝ effect of transposition along the antidiagonal.
       ∇¨                                 ⍝ For each element in the above vector, recurse.
      ∊                                   ⍝ Recursively flatten the results into a single
                                          ⍝ vector.


More details on "monadic transpose scan".

Dyalog documentation on error guards.