How do I compute the linear index of a 3D coordinate and vice versa?

Question

If I have a point (x, y z), how do I find the linear index, i for that point? My numbering scheme would be (0,0,0) is 0, (1, 0, 0) is 1, . . ., (0, 1, 0) is the max-x-dimension, .... Also, if I have a linear coordinate, i, how do I find (x, y, z)? I can't seem to find this on google, all the results are filled with other irrelevant stuff. Thank you!

Miaow

Are the coordinates always composed of integers? can you have negative coordinates? Do you have maximums for any axes besides the x axis? – Huguenot

Does each coordinate have the same number of divisions, or different? The last point is represented by (N,N,N) or (N1,N2,N3)? – Postal

Answer 1

There are a few ways to map a 3d coordinate to a single number. Here's one way.

some function f(x,y,z) gives the linear index of coordinate(x,y,z). It has some constants a,b,c,d which we want to derive so we can write a useful conversion function.

f(x,y,z) = a*x + b*y + c*z + d

You've specified that (0,0,0) maps to 0. So:

f(0,0,0) = a*0 + b*0 + c*0 + d = 0
d = 0
f(x,y,z) = a*x + b*y + c*z

That's d solved. You've specified that (1,0,0) maps to 1. So:

f(1,0,0) = a*1 + b*0 + c*0 = 1
a = 1
f(x,y,z) = x + b*y + c*z

That's a solved. Let's arbitrarily decide that the next highest number after (MAX_X, 0, 0) is (0,1,0).

f(MAX_X, 0, 0) = MAX_X
f(0, 1, 0) = 0 + b*1 + c*0 = MAX_X + 1
b = MAX_X + 1
f(x,y,z) = x + (MAX_X + 1)*y + c*z

That's b solved. Let's arbitrarily decide that the next highest number after (MAX_X, MAX_Y, 0) is (0,0,1).

f(MAX_X, MAX_Y, 0) = MAX_X + MAX_Y * (MAX_X + 1)
f(0,0,1) = 0 + (MAX_X + 1) * 0  + c*1 = MAX_X + MAX_Y * (MAX_X + 1) + 1
c = MAX_X + MAX_Y * (MAX_X + 1) + 1
c = (MAX_X + 1) + MAX_Y * (MAX_X + 1)
c = (MAX_X + 1) * (MAX_Y + 1)

now that we know a, b, c, and d, we can write your function as follows:

function linearIndexFromCoordinate(x,y,z, max_x, max_y){
    a = 1
    b = max_x + 1
    c = (max_x + 1) * (max_y + 1)
    d = 0
    return a*x + b*y + c*z + d
}

You can get the coordinate from the linear index by similar logic. I have a truly marvelous demonstration of this, which this page is too small to contain. So I'll skip the math lecture and just give you the final method.

function coordinateFromLinearIndex(idx, max_x, max_y){
    x =  idx % (max_x+1)
    idx /= (max_x+1)
    y = idx % (max_y+1)
    idx /= (max_y+1)
    z = idx
    return (x,y,z)
}

Huguenot

Great answer! I guess I'll just puzzle over your marvelous proof for 375 more years (but it makes sense now). Thanks a bunch. – Miaow

@Huguenot Hello! I realize that this question is almost 2 years old, but I was wondering: Would you maybe have a link to that math lecture you mentioned? Your method seems absolutely fantastic, so I was curious about the math behind it all. – Jermaine

you shouldn't alter code in edits after the answer is accepted - it shuold be done in a comment. – Zanze

@MarkAnderson, if I recall correctly, the math lecture only existed in my head while I was composing coordinateFromLinearIndex. I could edit the answer to show how I derived the final code block, but first I need to remember how I did it :-) – Huguenot

Can I interest you in the virtues of the divmod() function? I'm sure coordinateFromLinearIndex() could benefit from its unique properties. – Triliteral

Bookmarked! Great explanation of the math behind the algorithm. – Horwitz

@Huguenot is there any generic algorithm to make this conversion for higher dimensional arrays ? Let's say arrays with dimension 4 and 5 ? – Hawkes

For anyone who lands on this page seeking a "truly marvelous" demonstration, see my answer here: math.stackexchange.com/questions/3758576/…. – Baritone

Answer 2

If you have no upper limit on the coordinates, you can number them from origo and outwards. Layer by layer.

(0,0,0) -> 0
(0,0,1) -> 1
(0,1,0) -> 2
(1,0,0) -> 3
(0,0,2) -> 4
   :       :
(a,b,c) -> (a+b+c)·(a+b+c+1)·(a+b+c+2)/6 + (a+b)·(a+b+1)/2 + a

The inverse is harder, since you would have to solve a 3rd degree polynomial.

m1 = InverseTetrahedralNumber(n)
m2 = InverseTriangularNumber(n - Tetra(m1))
a = n - Tetra(m1) - Tri(m2)
b = m2 - a
c = m1 - m2

where

InverseTetrahedralNumber(n) = { x ∈ ℕ | Tetra(n) ≤ x < Tetra(n+1) } 
Tetra(n) = n·(n+1)·(n+2)/6 
InverseTriangularNumber(n) = { x ∈ ℕ | Tri(n) ≤ x < Tri(n+1) } 
Tri(n) = n·(n+1)/2

InverseTetrahedralNumber(n) could either be calculated from the large analytic solution, or searched for with some numeric method.

---

Here is my attempt at an algebraic solution (javascript). I am using the substitutions p = a+b+c, q = a+b, r = a to simplify the equations.

function index(a,b,c) {
    var r = a;
    var q = r + b;
    var p = q + c;
    return (p*(p+1)*(p+2) + 3*q*(q+1) + 6*r)/6;
}

function solve(n) {
    if (n <= 0) {
        return [0,0,0];
    }

    var sqrt = Math.sqrt;
    var cbrt = function (x) { return Math.pow(x,1.0/3); };

    var X = sqrt(729*n*n - 3);
    var Y = cbrt(81*n + 3*X);
    var p = Math.floor((Y*(Y-3)+3)/(Y*3));
    if ((p+1)*(p+2)*(p+3) <= n*6) p++;
    var pp = p*(p+1)*(p+2);

    var Z = sqrt(72*n+9-12*pp);
    var q = Math.floor((Z-3)/6);
    if (pp + (q+1)*(q+2)*3 <= n*6) q++;
    var qq = q*(q+1);

    var r = Math.floor((6*n-pp-3*qq)/6);
    if (pp + qq*3 + r*6 < n*6) r++;

    return [r, q - r, p - q];
}

Anglim