This answer complements Israel's answer (to which this edit was rejected).
The following expands his answer on how the L8 matrix was calculated. (I believe it is important, because the first result of googling with '2d Laplacian L8 matrix' directs to this link, however its calculation was missing.)
Previously, from the Taylor expansion
$$(1)\ \ f(x+\delta x,y+\delta y)=f(x,y)+(f_{x}\delta x + f_{y}\delta y)+ \frac{1}{2}(f_{xx}(\delta x)^{2}+2f_{xy} \delta x \delta y + f_{yy}(\delta y)^{2})+o((\delta x)^{3}+(\delta y)^{3})$$
One can get
$$(2) \qquad f(x+h,y)+f(x-h,y)=2f+f_{xx}h^{2}+o(h^{3})$$
$$(3) \qquad f(x,y+h)+f(x,y-h)=2f+f_{yy}h^{2}+o(h^{3})$$
By calculating (2)+(3):
$$(4) \ \ f(x,y+h)+f(x,y-h)+f(x+h,y)+f(x-h,y) = 4f +h^{2}
(f_{xx}+f_{yy})$$
Now, similar to (2) and (3), from the Tayler expansion (1) one can get:
$$(5)\ f(x+h,y+h) = f(x,y)+(f_x h+f_y h)+\frac{1}{2}(f_{xx}h^2+2f_{xy}h^2+f_{yy}h^2)+o(h^3)$$
$$(6)\ f(x+h,y-h) = f(x,y)+(f_x h-f_y h)+\frac{1}{2}(f_{xx}h^2-2f_{xy}h^2+f_{yy}h^2)+o(h^3)$$
$$(7)\ f(x-h,y+h) = f(x,y)+(-f_x h+f_y h)+\frac{1}{2}(f_{xx}h^2-2f_{xy}h^2+f_{yy}h^2)+o(h^3)$$
$$(8)\ f(x-h,y-h) = f(x,y)+(-f_x h-f_y h)+\frac{1}{2}(f_{xx}h^2-2f_{xy}h^2+f_{yy}h^2)+o(h^3)$$
By summing up (5),(6),(7),(8) one can get:
$$(9)\ \ f(x+h,y+h)+f(x+h,y-h)+f(x-h,y+h)+f(x-h,y-h)=4f+2h^2(f_{xx}+f_{yy})$$
By adding (4) and (9) together,
$$(10)\ \ f(x,y+h)+f(x,y-h)+f(x+h,y)+f(x-h,y)+f(x+h,y+h)+f(x+h,y-h)+f(x-h,y+h)+f(x-h,y-h)=8f+3h^2(f_{xx}+f_{yy})=8f+3h^2\nabla^{2}f$$
Therefore,
$$\nabla^{2}f = \frac{1}{3h^2}(-8f+f(x,y+h)+f(x,y-h)+f(x+h,y)+f(x-h,y)+f(x+h,y+h)+f(x+h,y-h)+f(x-h,y+h)+f(x-h,y-h))$$
Letting $h=1$, The matrix form becomes$$
L8= \frac{1}{3h^2}\begin{pmatrix}
c_{f(x-h,y+h)} & c_{f(x,y+h)} & c_{f(x+h,y+h)}\\
c_{f(x-h,y)} & c_{f(x,y)} & c_{f(x+h,y)}\\
c_{f(x-h,y-h)} & c_{f(x,y-h)} & c_{f(x+h,y-h)}\\
\end{pmatrix} = \frac{1}{3}\begin{pmatrix}
1 & 1 & 1\\
1 & -8 & 1\\
1 & 1 & 1\\
\end{pmatrix}
$$
In addition, one often see a more accurate kernel from the literatures, which can be calculated as the following:
4x(4)+(9):
$$(11)\ \ 4f(x,y+h)+4f(x,y-h)+4f(x+h,y)+4f(x-h,y)+f(x+h,y+h)+f(x+h,y-h)+f(x-h,y+h)+f(x-h,y-h)=20f+6h^2\nabla^{2}f$$
Therefore,
$$\nabla^{2}f = \frac{1}{6h^2}(-20f+4f(x,y+h)+4f(x,y-h)+4f(x+h,y)+4f(x-h,y)+f(x+h,y+h)+f(x+h,y-h)+f(x-h,y+h)+f(x-h,y-h))$$
Letting $h=1$, The matrix form becomes$$
\hat{L}8= \frac{1}{6h^2}\begin{pmatrix}
c_{f(x-h,y+h)} & c_{f(x,y+h)} & c_{f(x+h,y+h)}\\
c_{f(x-h,y)} & c_{f(x,y)} & c_{f(x+h,y)}\\
c_{f(x-h,y-h)} & c_{f(x,y-h)} & c_{f(x+h,y-h)}\\
\end{pmatrix} = \frac{1}{6}\begin{pmatrix}
1 & 4 & 1\\
4 & -20 & 4\\
1 & 4 & 1\\
\end{pmatrix}
$$