In the paper about Speeded-Up Robust Features, it says that in order to localize points, interpolation of nearby data is needed to find the location in space and scale. This is done by fitting a 3D quadratic. In order to do this, the determinant of the Hessian function H(x, y, sigma) is expressed as a Taylor expansion up to quadratic terms centered at detected location.

My question is, how do I solve for the elements in the matrices shown? How do I solve for dx, dxx, dxy, etc.? Thank you!