Why does Shamir secret sharing appear to need ordered shares?

The implementation of Shamir secret sharing in this code, only generates the original image if the shares are provided in consecutive order (ex: [2,3,4]) and won't work in any other share order (ex: [2,4,6] or [4,1,3]). However, Shamir secret reconstruction does not require the shares to be in any order, then why does this fail?

import numpy as np
from scipy.interpolate import lagrange as lag
import os

def polynomial(img, n, k):
    '''
        Generate lagrange polynomial of degree k-1
        f(x) = c1(x^k) + c2(x^k-1) .... + secret (mod prime)
        prime > secret
    '''
    coef = np.random.randint(low = 0, high = 251, size = (img.shape[0], k - 1)) #Coefficients should not exceed value of prime number chosen
    gen_imgs = []
    for i in range(1, n + 1):
        base = np.array([i ** j for j in range(1, k)])
        base = np.matmul(coef, base)
        imgValue_ = (img + base) % 251
        gen_imgs.append(imgValue_)
    return np.array(gen_imgs)

def reconstruct(imgs, index, k):
    '''
        Reconstruct image using share index values for k shares
    '''
    print("Shares: ", index)
    assert imgs.shape[0] >= k
    x = np.array(index)
    dim = imgs.shape[1]
    img = []
    for i in range(dim):
        if i % 10000 == 0:
            print("Reconstructing pixel ", i, " of ", dim, " pixels")
        y = imgs[:, i]
        poly = lag(x, y)
        pixel = poly(0) % 251
        img.append(pixel)
    return np.array(img)

    
if __name__ == "__main__":
    pathPrefix = path.split('.')[0]
    os.makedirs(pathPrefix, exist_ok=True)

    img_flattened, shape = util.read_image(path)
    gen_imgs = polynomial(img_flattened, n = n, k = k)
    to_save = gen_imgs.reshape(n, *shape)
    for i, img in enumerate(to_save):
        Image.fromarray(img.astype(np.uint8)).save(pathPrefix + "/share" + "_{}.jpeg".format(i + 1))
    
    #Secret reconstruction
    shareIndex = list(map(int, input("\nEnter the index of shares to reconstruct: ").split()))
    origin_img = reconstruct(gen_imgs, shareIndex, k = k)
    origin_img = origin_img.reshape(*shape)
    Image.fromarray(origin_img.astype(np.uint8)).save(pathPrefix + "/reconstructed_image.jpeg")

The read_image function has the following definition

def read_image(path):
    '''
        Reads image from file and converts it into numpy array in greyscale
    '''
    img = Image.open(path).convert('L') 
    img_array = np.asarray(img)
    return img_array.flatten(), img_array.shape

The implementation of lag is

def lagrange(x, w):
    M = len(x)
    p = poly1d(0.0)
    for j in range(M):
        pt = poly1d(w[j])
        for k in range(M):
            if k == j:
                continue
            fac = x[j]-x[k]
            pt *= poly1d([1.0, -x[k]])/fac
        p += pt
    return p

However, Shamir secret reconstruction does not require the shares to be in any order, then why does this fail?

Well, Shamir's secret sharing works in a finite field; in your code, you use the field $Z/251$ (which is a valid field).

However, in your reconstruction code, you have the line:

        pt *= poly1d([1.0, -x[k]])/fac

To be correct, the division here needs to be a division in the field $Z/251$ [1]; instead, what this line will do is division over the reals.

Now, the obvious question is: why does it work in the $[2,3,4]$ case? Well, it might be mostly a coincidence; in the $3-2$ and the $4-3$ case, dividing by 1 will give the correct result. And, in the $4-2$ case, that will happen to give the correct result if poly1d([1.0, -x[k]]) happens to be even (which, if that's a random value, will happen half the time).

[1]: In $Z/251$, the value $a/b$ is that number $c$ such that $b*c \equiv a \pmod {251}$; if $b \ne 0$, it will always exist (and be unique). To compute it, well, the easiest way in this case probably would be to construct a table of values $1/b$ for the 250 possible values of $b$, and use the identity $a/b = a * (1/b)$ (not constant time, however we probably don't care about constant time in this case)

The shares are really pairs $(x_i,f(x_i))$ where $x_i$ is the coordinate and $f(x_i)$ is its evaluation under the lagrange interpolation polynomial. Most likely your lists are $[f(x_i): 1\leq i\leq 3]$ with some standard point coordinates $[x_i:1\leq i\leq 3]$ assumed.

In that sense, order does matter. Pixels in an image are ordered so why wouldn't it matter?

Also note that in the code fragment, no details of the reconstruction algorithm (lagrange interpolation) is given.

Edit: The pixel location in the image is fixed. The obvious way of doing this would have that location (coordinate if you will) encode the value $x_i,$ while the scalar representing the level of the pixel (let's pretend it's in grayscale) as an integer would represent the value $f(x_i)$. As obvious from the code the points $x_i$ are from the array (ordered list) range(M) presumably $[0,1,\ldots,M-1]$ (if I've interpreted the command "range" correctly, maybe it is $[1,\ldots,M].$ So the $x_i$ are fixed and ordered.