Suppose f is a function in x and y. The bi-Laplacian of f is
D (D f) = fxxxx + 2fxxyy + fyyyy where fxxxx is the fourth derivative of x wrt x.

We want to express the bi-Laplacian of f at a point x as a linear combination of
the value of f at nearby points. In 2D we can use the following 13-point stencil.

 2        *
 1     *  *  *
 0  *  *  x  *  *
-1     *  *  *
-2        *
   -2 -1  0  1  2

The point in question is marked as x at the center (0, 0).
fxxxx(0, 0) can be found by repeatedly applying the central difference formula
for the second derivative in x:

fxx(-1, 0) = f(0, 0) + f(-2, 0) - 2f(-1, 0)
fxx( 0, 0) = f(1, 0) + f(-1, 0) - 2f( 0, 0)
fxx( 1, 0) = f(0, 0) + f( 2, 0) - 2f( 1, 0)

so
fxxxx(0, 0) = fxx(-1, 0) + fxx(1, 0) - 2fxx(0, 0)
            = f(-2, 0) - 4f(-1, 0) + 6f(0, 0) - 4f(1, 0) + f(2, 0) (1)
This corresponds to a 5-point mask:
 0  1 -4  6 -4  1
   -2 -1  0  1  2
            
We get a similar formula for fyyyy(0, 0):
fyyyy(0, 0) = f(0, -2) - 4f(0, -1) + 6f(0, 0) - 4f(0, 1) + f(0, 2) (2)
The mask for this one is:
 2  1
 1 -4
 0  6
-1 -4
-2  1
    0

To find a formula for fxxyy(0, 0), we apply the central difference formula for
the second derivative in x first, then y:

fxx(0,  1) = f(-1,  1) + f(1,  1) - 2f(0,  1)
fxx(0,  0) = f(-1,  0) + f(1,  0) - 2f(0,  0)
fxx(0, -1) = f(-1, -1) + f(1, -1) - 2f(0, -1)

so
fxxyy(0, 0) = fxx(0, 1) + fxx(0, -1) - 2fxx(0, 0)
            =    f(-1, -1) - 2f(0, -1) +  f(1, -1)
              - 2f(-1,  0) + 4f(0,  0) - 2f(1,  0) 
              +  f(-1,  1) - 2f(0,  1) +  f(1,  1) (3)
This formula gives us the following 9-point mask (the left column and the bottom
row are the indices and not part of the mask):
 1  1 -2  1
 0 -2  4 -2
-1  1 -2  1
   -1  0  1
              
Plugging in the formulas (1), (2), (3) into D (D f) = fxxxx + 2fxxyy + fyyyy
we get the following mask for the bi-Laplacian at (0, 0):

 2        1
 1     2 -8  2
 0  1 -8 20 -8  1
-1     2 -8  2
-2        1
   -2 -1  0  1  2
   
By simply applying the same mask in each plane xy, yz, zx we can generalize this
scheme to 3D, where the bi-Laplacian of f is
fxxxx + fyyyy + fzzzz + 2fxxyy + 2fyyzz + 2fzzxx

At the boundaries, we notice that the masks derived in (1), (2) and (3) each can
be slided back into the grid independently in x and y should they get outside
the grid. For example, the mask in (3) can be translated up and to the right for
a point near the lower left corner of the grid as follows:

 2  1 -2  1
 1 -2  4 -2
 0  1 -2  1
    0  1  2
    
We still have a 9-point mask but the point in question (0, 0) is now at the
lower left of the mask. The masks generated by (1) and (2) must also be
translated accordingly, and the summation of all three masks looks like this:

 4  *
 3  *
 2  *  *  *
 1  *  *  *
 0  x  *  *  *  *
    0  1  2  3  4
    
In code, it would probably be easier to keep the masks for (1), (2) and (3)
separated, translate them independently and add them up at the end.