kernel to calculate second order partial derivative of digital image

im working on image stacks, and I need to calculate second order partial derivatiatives of it. I already know how to calculate derivative on x and y axis, using coefficients calculator: https://web.media.mit.edu/~crtaylor/calculator.html

So for example to calculate five point second order derivative on x axis we have formula like:

So I already know how to calculate derivatives on x and y axes, but I also need to get formulas how to calculate first and second order partial derivatives for f'xy and f''xy. I know that I could simply calculate derivatives on x axis first, then calculate from it on y axis, but I would like to obtain simplified formula that uses x and y pixels in image to calculate f'xy/f''xy. Could anyone explain me how to obtain them?

Update: I found somewhere on the web the formula for second order f''xy derivative for images that looks like that: (inX[2]-1inX[0]+inY[2]-1inY[0])/4;
Would like also to get an explanation whether its correct, and how its acquired

Your five-point derivative kernel is a 1D kernel. Applied along the x axis it gives the partial derivative for x, applied along the y axis it gives the partial derivative for y.

The $$\frac{\partial^2}{\partial x \partial y}$$ derivative would need a 2D square kernel. It is more efficient to compute this by applying two first order partial derivatives.

$$\frac{\partial^2}{\partial x \partial y} f = \left( \frac{\partial}{\partial x} \frac{\partial}{\partial x} \right) f = \frac{\partial}{\partial x} \left( \frac{\partial}{\partial x} f \right)$$

Since $$\frac{\partial}{\partial x} f$$ is a convolution,

$$\frac{\partial}{\partial x} f = k_x \ast f$$

we can write the earlier equation as

$$\frac{\partial^2}{\partial x \partial y} f = \left( k_x \ast k_y \right) \ast f = k_x \ast \left( k_y \ast f \right)$$

So we can say that the 2D square convolution kernel that you are looking for is $$k_x \ast k_y$$ (with $$k_x$$ the kernel that computes the x-derivative and $$k_y$$ the one that computes the y-derivative). Note that, if you think of your 1D kernels as vectors, this convolution is equivalent to the outer product of the two vectors.

However, it is more expensive to compute a convolution with one 5x5 kernel (25 multiplication and additions per pixel) than to compute two convolutions with 5x1 and 1x5 kernels (10 multiplications and additions per pixel). I would compute two first order derivatives in sequence.