How was the 2D discrete Laplacian matrix calculated? If I convolute a 2D signal with the L4 kernel, I get the discrete Laplacian. But how was the L4 kernel calculated?
Wikipedia isn't very clear on this and I haven't found other information.
(https://en.wikipedia.org/wiki/Discrete_Laplace_operator#Implementation_via_operator_discretization)
I've also read this but it's still unclear how the convolution of [1, -2, 1] is L4 transposed.
L4 being this:
 L4 =
  [ 0 -1  0;
   -1  4 -1;
    0 -1  0];

EDIT:
Also, how is the L8 produced?
 L8 =
  [-1 -1 -1;
   -1  8 -1;
   -1 -1 -1];

 A: Basically the 2D "L4" discrete Laplacian operator is constructed by using 4 surrounding points from a central stencil point. These stencil points are north, south, east and west from the central point.
In space these 5 stencil points are then given by:
$$ (x,y), (x+h,y), (x,y+h), (x-h,y), (x,y-h)$$
Remembering that
$$\nabla^{2} f = f_{xx}+f_{yy} $$
and that a 2D Taylor expansion (up to a point) may be seen as:
$$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})$$
We can then find that
$$f(x\pm h, y)= f(x,y) \pm f_{x}h+\frac{1}{2}f_{xx}h^{2}+O(h^{3})$$
$$f(x, y\pm h)= f(x,y) \pm f_{y}h+\frac{1}{2}f_{yy}h^{2}+O(h^{3})$$
By summing expressions we may obtain
$$(1) \qquad f(x+h,y)+f(x-h,y)=2f+f_{xx}h^{2}+O(h^{3})$$
$$(2) \qquad f(x,y+h)+f(x,y-h)=2f+f_{yy}h^{2}+O(h^{3})$$
Calculating (1) + (2):
$$ f(x,y+h)+f(x,y-h)+f(x+h,y)+f(x-h,y) = 4f +h^{2}(f_{xx}+f_{yy})$$
Then using our very first equation (Cartesian Laplacian definition) we are certain that we approximated our Laplacian up to 2nd order:
$$\nabla^{2} f = f_{xx}+f_{yy}= \frac{1}{h^{2}} \Big[f(x,y+h)+f(x,y-h)+f(x+h,y)+f(x-h,y) -4f(x,y) \Big]$$
As indeed we want to see this in a discrete 2D space, we may correspond the coefficients of our stencils onto a matrix where $h=1$
$$
    L4 = \frac{1}{h^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} = \begin{pmatrix}
    0 & 1 & 0\\
    1 & -4 & 1\\
    0 & 1 & 0\\
    \end{pmatrix}
$$
One can follow a similar procedure to obtain L8 by using a 9-stencil configuration where the stencils would be
$$  f(x,y), f(x+h,y+h),f(x+h,y),f(x+h,y-h), f(x,y-h), f(x,y+h), f(x-h,y-h),f(x-h,y+h),f(x-h,y) $$
A: 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}
$$
