3d array traversal originating from center

Question

I'm trying to find a traversal order for a 3d array with uniform dimension n.The traversal order should hereby be sorted ascending by it's distance to the cube's center (order of cells with equal indices is arbitrary). Example for a 2d-Array:

7 4 8
3 0 1
6 2 5

Which basically is the distance in Manhattan metric:

2 1 2
1 0 1
2 1 2

Traversal as relative coordinates with respect to the origin:

[ 0, 0]
[ 1, 0]
[ 0,-1]
[-1, 0]
[ 0, 1]
[ 1,-1]
[-1,-1]
[-1, 1]
[ 1, 1]

I know of some ways to solve this, for instance precalculating all indices and sorting them according to their distance to the origin. However, since this algorithm is intended to perform on GPU, there are some restrictions:

1. No recursion (I'm aware of the possibility resolving recursion into an iterative algorithm - maintaining a stack is however not a suitable solution in my experience)
2. No offline calculation (= calculation on the CPU and transferring the result to the GPU). The solution needs to be as flexible as possible

While searching for solutions I stumbled upon this question, which is exactly the problem I tend to solve, the accepted answer albeit envolves a tree structure, which does not fit the specified requirements: 3D Array Traversal in Different Order

I also thought of a way to create the indices using spherical coordinates, which unfortunately does not yield the correct order. What is an appropriate algorithm to generate the given traversal order for 3d-arrays?

Edit: Stormwind provided an excellent alternative description for the given Problem: "[The problem] is actually about converting addressing from one space to another. Converting between a 1-dimensional and 2-dimensional is simple, like 1,2,3,... to (1,1),(1,2)(2,1)... But this is more like converting from an ascending 1-dimensional (or at least square) to an "ascending octahedron layered" space, when "ascending" means "innermost layer first" in addition to an existing (though arbitrary) incrementing order at each layer surface."

Answer 1

After thinking for a while I've came up with an idea to represent 3d array as a sequence of nodes with directions: +i, -i, +j, -j, +k, -k.

Approach

For 2-dimensional array it would be sufficient to have only three rules:

1. Each iteration over each node moves it along its axis in its direction, i.e. node +j will increment 2nd index, node -i will decrement 1st index.
2. There are two kinds of nodes: Main and Secondary. Main axis have one index 0. Each iteration over Main i and j axis nodes (I'll call them I and J) produces Secondary node rotated clockwise 90 degrees:
   • +I -> -j
   • -J -> -i
   • -I -> +j
   • +J -> +i
3. Each node has lifetime which decrements every iteration. Lifetime of the node equals (n-1)/2 for odd values of n (for even values see code below). After lifetime goes to 0 the node should be removed.

To enable 3rd dimension another rule should be applied:

4. Third type of nodes with direction along k axis (here - depth) produces set of I and J axis in each iteration:
   • +K -> +I, -J, -I, +J
   • -K -> +I, -J, -I, +J

Here is how it will look:

With such approach elemnts will be automatically sorted by Manhattan distance as in Arturo Menchaca solution.

Implementation

Here is python code that does what I've described. There is a lot of space for improvements this is just a proof of concept. It has no sorting, no recursion and I don't see any offline calculations. It also contains few tests. Run

NO = ( 0, 0, 0, 2, 0)
Pi = (+1, 0, 0, 0, 0)
PI = (+1, 0, 0, 0, 1)
Pj = ( 0,+1, 0, 0, 0)
PJ = ( 0,+1, 0, 0, 1)
PK = ( 0, 0,+1, 0, 2)
Mi = (-1, 0, 0, 1, 0)
MI = (-1, 0, 0, 1, 1)
Mj = ( 0,-1, 0, 1, 0)
MJ = ( 0,-1, 0, 1, 1)
MK = ( 0, 0,-1, 1, 2)
#      i  j  k  ^  ^
#               |  Id for comparison
#               Lifetime index

PRODUCE = {
    PI: [ Mj ], # +I -> -j
    MJ: [ Mi ], # -J -> -i
    MI: [ Pj ], # -I -> +j
    PJ: [ Pi ], # +J -> +i
    NO: [ NO ],
    Pi: [ NO ],
    Pj: [ NO ],
    Mi: [ NO ],
    Mj: [ NO ],
    PK: [ PI, MI, PJ, MJ ], # +K -> +I, -J, -I, +J
    MK: [ PI, MI, PJ, MJ ], # -K -> +I, -J, -I, +J
}


class Index:
    LastDistance = 0
    NumberOfVisits = 0
    MinIndex = 0
    MaxIndex = 0
    def __init__(self, i, j, k, lifetime, direction):
        self.globalLifetime = lifetime
        self.direction = direction

        # Assign parent's position
        self.i = i
        self.j = j
        self.k = k

        # Step away from parent
        self.lifetime = lifetime[direction[3]]
        self.step()

    def isLive(self):
        return self.lifetime > 0

    def visit(self):
        Index.NumberOfVisits += 1
        distance = self.distance()
        if distance < Index.LastDistance:
           raise NameError("Order is not preserved")
        Index.LastDistance = distance
        Index.MinIndex = min(self.i, Index.MinIndex)
        Index.MinIndex = min(self.j, Index.MinIndex)
        Index.MinIndex = min(self.k, Index.MinIndex)
        Index.MaxIndex = max(self.i, Index.MaxIndex)
        Index.MaxIndex = max(self.j, Index.MaxIndex)
        Index.MaxIndex = max(self.k, Index.MaxIndex)
        print("[{}, {}, {}]".format(self.i, self.j, self.k))

    def step(self):
        # Move in your direction
        self.i += self.direction[0]
        self.j += self.direction[1]
        self.k += self.direction[2]

    def iterate(self):
        self.lifetime -= 1

    def produce(self, result):
        for direction in PRODUCE[self.direction]:
            self.create(direction, result)

    def create(self, direction, result):
        index = Index(self.i, self.j, self.k, self.globalLifetime, direction)
        if index.isLive():
            result.append(index)

    def distance(self):
        # Manhattan Distance
        return abs(self.i) + abs(self.j) + abs(self.k)

def Traverse(N):
    TotalNumber = N*N*N
    halfA = halfB = (N-1)/2
    if N % 2 == 0:
        halfA = N/2
        halfB = N/2-1

    MinIndex = -min(halfB, halfA)
    MaxIndex = max(halfB, halfA)

    lifetime = (halfA, halfB, 0)

    SecondaryNodes = []
    MainNodes = []
    KNodes = []

    # visit center
    center = Index(0, 0, 0, lifetime, NO)
    center.visit()

    center.create(PI, MainNodes)
    center.create(MI, MainNodes)
    center.create(PJ, MainNodes)
    center.create(MJ, MainNodes)
    center.create(PK, KNodes)
    center.create(MK, KNodes)

    while len(SecondaryNodes) + len(MainNodes) + len(KNodes) > 0:

        # First - visit all side nodes
        temp = []
        for m in SecondaryNodes:
            m.visit()
            m.step()
            m.iterate()
            # Save node only if it is alive
            if m.isLive():
                temp.append(m)

        SecondaryNodes = temp

        # Second - visit all main nodes as they may produce secondary nodes
        temp = []
        for m in MainNodes:
            m.visit() # 1 - Visit
            m.produce(SecondaryNodes) # 2 - Produce second
            m.step() # 3 - Step 
            m.iterate() # 4 - Lose a life
            if m.isLive():
                temp.append(m)

        MainNodes = temp

        # Third - visit all K nodes as they may produce main nodes
        temp = []
        for m in KNodes:
            m.visit()
            m.produce(MainNodes)
            m.step()
            m.iterate()
            if m.isLive():
                temp.append(m)

        KNodes = temp
    if TotalNumber != Index.NumberOfVisits:
        raise NameError("Number of visited elements does not match {}/{}".format(Index.NumberOfVisits, TotalNumber))
    if MinIndex != Index.MinIndex:
        raise NameError("Minimal index is out of bounds {}/{}".format(Index.MinIndex, MinIndex))
    if MaxIndex != Index.MaxIndex:
        raise NameError("Maximal index is out of bounds {}/{}".format(Index.MaxIndex, MaxIndex))

Traverse(6)

Implementation Simplified

Helper class to store index:

class Index:
    def __init__(self, i, j, k, lifetime):
        self.i = i
        self.j = j
        self.k = k
        self.lifetime = lifetime

    def visit(self):
        print("[{}, {}, {}]".format(self.i, self.j, self.k))

Set of functions to iterate Main nodes in proper direction:

def StepMainPlusI(mainPlusI, minusJ, lifetime):
    result = []
    for m in mainPlusI:
        if lifetime > 0:
            minusJ.append(Index(m.i, m.j-1, m.k, lifetime))
        m.lifetime -= 1
        m.i += 1
        if m.lifetime > 0:
            result.append(m)
    return result

def StepMainMinusJ(mainMinusJ, minusI, lifetime):
    result = []
    for m in mainMinusJ:
        if lifetime > 0:
            minusI.append(Index(m.i-1, m.j, m.k, lifetime))
        m.lifetime -= 1
        m.j -= 1
        if m.lifetime > 0:
            result.append(m)
    return result

def StepMainMinusI(mainMinusI, plusJ, lifetime):
    result = []
    for m in mainMinusI:
        if lifetime > 0:
            plusJ.append(Index(m.i, m.j+1, m.k, lifetime))
        m.lifetime -= 1
        m.i -= 1
        if m.lifetime > 0:
            result.append(m)
    return result

def StepMainPlusJ(mainPlusJ, plusI, lifetime):
    result = []
    for m in mainPlusJ:
        if lifetime > 0:
            plusI.append(Index(m.i+1, m.j, m.k, lifetime))
        m.lifetime -= 1
        m.j += 1
        if m.lifetime > 0:
            result.append(m)
    return result

Set of functions to iterate a third dimensional K nodes:

def StepMainPlusK(mainPlusK, mainPlusI, mainMinusI, mainPlusJ, mainMinusJ, lifetimeA, lifetimeB):
    result = []
    for m in mainPlusK:
        if lifetimeA > 0:
            mainPlusI.append(Index(+1, 0, m.k, lifetimeA))
            mainPlusJ.append(Index(0, +1, m.k, lifetimeA))
        if lifetimeB > 0:
            mainMinusI.append(Index(-1, 0, m.k, lifetimeB))
            mainMinusJ.append(Index(0, -1, m.k, lifetimeB))
        m.lifetime -= 1
        m.k += 1
        if m.lifetime > 0:
            result.append(m)
    return result

def StepMainMinusK(mainMinusK, mainPlusI, mainMinusI, mainPlusJ, mainMinusJ, lifetimeA, lifetimeB):
    result = []
    for m in mainMinusK:
        if lifetimeA > 0:
            mainPlusI.append(Index(+1, 0, m.k, lifetimeA))
            mainPlusJ.append(Index(0, +1, m.k, lifetimeA))
        if lifetimeB > 0:
            mainMinusI.append(Index(-1, 0, m.k, lifetimeB))
            mainMinusJ.append(Index(0, -1, m.k, lifetimeB))
        m.lifetime -= 1
        m.k -= 1
        if m.lifetime > 0:
            result.append(m)
    return result

These two functions have two different lifetime parameters for the case when n is odd and one half can be less than another. I've divided them by sign - negatively oriented will have lower half of indexes.

Set of functions to iterate Secondary nodes:

def StepPlusI(plusI):
    result = []
    for m in plusI:
        m.i += 1
        m.lifetime -= 1
        if m.lifetime > 0:
            result.append(m)
    return result

def StepMinusI(minusI):
    result = []
    for m in minusI:
        m.i -= 1
        m.lifetime -= 1
        if m.lifetime > 0:
            result.append(m)
    return result

def StepPlusJ(plusJ):
    result = []
    for m in plusJ:
        m.j += 1
        m.lifetime -= 1
        if m.lifetime > 0:
            result.append(m)
    return result

def StepMinusJ(minusJ):
    result = []
    for m in minusJ:
        m.j -= 1
        m.lifetime -= 1
        if m.lifetime > 0:
            result.append(m)
    return result

And the main function:

def Traverse(N):
    halfA = halfB = (N-1)/2
    if N % 2 == 0: # size is even
        halfA = N/2
        halfB = N/2-1

    # visit center
    Index(0,0,0,0).visit()

    # Secondary nodes
    PlusI  = []
    MinusI = []
    PlusJ  = []
    MinusJ = []

    # Main nodes
    MainPlusI  = []
    MainMinusI = []
    MainPlusJ  = []
    MainMinusJ = []
    MainPlusK  = []
    MainMinusK = []

    # Add Main nodes
    if halfA > 0:
        MainPlusI.append(  Index(+1, 0, 0, halfA) )
        MainPlusJ.append(  Index(0, +1, 0, halfA) )
        MainPlusK.append(  Index(0, 0, +1, halfA) )

    if halfB > 0:
        MainMinusI.append( Index(-1, 0, 0, halfB) )
        MainMinusJ.append( Index(0, -1, 0, halfB) )
        MainMinusK.append( Index(0, 0, -1, halfB) )

    # Finish condition flag
    visited = True
    while visited:
        visited = False

        # visit all Main nodes
        for m in MainPlusI:
            m.visit()
            visited = True
        for m in MainMinusI:
            m.visit()
            visited = True
        for m in MainPlusJ:
            m.visit()
            visited = True
        for m in MainMinusJ:
            m.visit()
            visited = True
        for m in MainPlusK:
            m.visit()
            visited = True
        for m in MainMinusK:
            m.visit()
            visited = True

        # Visit all Secondary nodes
        for m in PlusI:
            m.visit()
            visited = True
        for m in MinusI:
            m.visit()
            visited = True
        for m in PlusJ:
            m.visit()
            visited = True
        for m in MinusJ:
            m.visit()
            visited = True

        # Iterate Secondary nodes first
        PlusI = StepPlusI(PlusI)
        MinusI = StepMinusI(MinusI)
        PlusJ = StepPlusJ(PlusJ)
        MinusJ = StepMinusJ(MinusJ)

        # Iterate all Main nodes as they might generate Secondary nodes
        MainPlusI = StepMainPlusI(MainPlusI, MinusJ, halfB)
        MainMinusJ = StepMainMinusJ(MainMinusJ, MinusI, halfB)
        MainMinusI = StepMainMinusI(MainMinusI, PlusJ, halfA)
        MainPlusJ = StepMainPlusJ(MainPlusJ, PlusI, halfA)

        # Iterate K nodes last as they might produce Main nodes
        MainPlusK = StepMainPlusK(MainPlusK, MainPlusI, MainMinusI, MainPlusJ, MainMinusJ, halfA, halfB)
        MainMinusK = StepMainMinusK(MainMinusK, MainPlusI, MainMinusI, MainPlusJ, MainMinusJ, halfA, halfB)

And the live example Code