Implementing an L-System conversion from R2 to R1 for symmetric, square matrices

Question

I have been reading for some time about L-Systems, and specifically the Hilbert Space filling curve. I am interested in writing a function to convert upper-triangular matrix coordinates into an ordered, 1-dimensional space filling curve.

The usual Hilbert space filling curve is derived with a single, first-order curve:

Then

I have successfully implemented the mapping algorithm listed on Wikipedia,

//rotate/flip a quadrant appropriately
void rot(int n, int *x, int *y, int rx, int ry) {
    if (ry == 0) {
        if (rx == 1) {
            *x = n-1 - *x;
            *y = n-1 - *y;
        }

        //Swap x and y
        int t  = *x;
        *x = *y;
        *y = t;
    }
}

//convert (x,y) to d
int xy2d (int n, int x, int y) {
    int rx, ry, s, d=0;
    for (s=n/2; s>0; s/=2) {
        rx = (x & s) > 0;
        ry = (y & s) > 0;
        d += s * s * ((3 * rx) ^ ry);
        rot(n, &x, &y, rx, ry);
    }
    return d;
}

//convert d to (x,y)
void d2xy(int n, int d, int *x, int *y) {
    int rx, ry, s, t=d;
    *x = *y = 0;
    for (s=1; s<n; s*=2) {
        rx = 1 & (t/2);
        ry = 1 & (t ^ rx);
        rot(s, x, y, rx, ry);
        *x += s * rx;
        *y += s * ry;
        t /= 4;
    }
}

However, I need to add two more first order curves to this process: one that is a right angle and the other which is a wedge. These two new curves do not need to be rotated or flipped, and in combination with the usual Hilbert 1st order curve, one can fill an upper triangular matrix (see my sketch below):

That is, the upper left of the red right-angle again becomes a right angle, the upper right becomes a Hilbert curve, and the lower right part of the right-angle becomes a wedge, and so on.

I would like to implement the above proposed L-System conversion, but every single tutorial on L-Systems is about STARTING AT ZERO and then DRAWING A LINE that fills the space without crossing itself (i.e. writing a string of steps to trace the whole space, instead of converting a subset of coordinates in the space to there 1-D counterpart).

Can anyone offer any intuition on how I might incorporate these new restriction into the existing code? Are there any resources that explain L-System conversions (as opposed to drawing systems), and how to construct them in code??

Answer 1

I've defined the following function in R

State 0 indicates that one is located within a right-angle region

State 1 indicates that one is located within a Hilbert curve

State 2 indicates that one is located within a wedge region

rot <- function(n, x, y, rx, ry) {
  if (ry == 1) {
    if (rx == 0) {
      x <<- n-1 - x
      y <<- n-1 - y
    }

    temp <- x
    x <<- y
    y <<- temp
  }
}

xy2d <- function(n, x, y) {
  # rx = x region (right = 1, left equals 0) 
  # ry = y region (top = 1, bottom = 0)
  # s = half the dimension
  # d = linear coordinate value
  d=0L
  state = 0L
  for (s in 2L^(log2(n/2L):0)) {
    rx = bitwAnd(x, s) > 0
    ry = bitwAnd(y, s) > 0
    if (state == 0) {
      # *--2  | 0  ->  *--2   | 1 -> 0--3   | 2 -> *--2   |  
      # 0--1  |        0--1   |      1--2   |      1--0   |

      # if x on right region, add lower left triangle (region 0)
      #      (s * (s + 1) / 2) * rx
      # if y on upper region, add lower right square as well (both regions 0 and 1)
      #     (s * s) * ry
      d = d + (s * (s + 1) / 2) * rx + (s * s) * ry
      d
      # change the state
      #  rx  ry | (rx + ry)
      #  0   0  |  0
      #  1   0  |  1
      #  1   1  |  2
      state = rx + ry
      state
    }else if (state == 1) {
      # 0--3  | 0  ->  0--1   | 1,2 -> 0--3   | 3 -> 2--3   |  
      # 1--2  |        3--2   |        1--2   |      1--0   |


      # Determine how many squares to add
      # equals number in s-by-s square times it's order in the curve
      #  rx  ry | (3 * rx) ^ (1 - ry)
      #  0   1  |  0 ^ 0 = [00] ^ [00] = 00 = 0
      #  0   0  |  0 ^ 1 = [00] ^ [01] = 01 = 1
      #  1   0  |  3 ^ 1 = [11] ^ [01] = 10 = 2
      #  1   1  |  3 ^ 0 = [11] ^ [00] = 11 = 3
      d = d + s * s * bitwXor(3 * rx, 1 - ry)
      d
      # ROTATE
      rot(n, x, y, rx, ry)
      x
      y

    }else if (state == 2) {
      # *--2  | 0  ->  2--1   | 1 -> *--2   | 3 -> *--2   |  
      # 1--0  |        3--0   |      1--0   |      0--1   |

      d = d + (s * s) * (rx == ry) + (s * (s + 1) / 2) * (rx == 1 && ry == 1)

      #  rx  ry |  Target | state - (rx + ry)
      #_________|_________|__________________
      #  1   0  |  1      | 2 - (1 + 0) = 1
      #  0   0  |  2      | 2 - (0 + 0) = 2
      #  1   1  |  0      | 2 - (1 + 1) = 0

      # new_state = current_state - (rx + ry)
      state = 2 - (rx + ry)
      if (state == 1) {
        # If the state goes from state-2 to state-1, then rotate the hilbert 
        # curve 180 degrees (i.e. reflect horizonally and vertically)
        x <- n-1 - x
        y <- n-1 - y
      }
    }
  }
  return(d)
}

xy2d(8,5,0)

[1] 16
```