Python - Find x and y values of a 2D gaussian given a value for the function

Question

I have a 2D gaussian function f(x,y). I know the values x₀ and y₀ at which the peak g₀ of the function occurs. But then I want to find xₑ and yₑ values at which f(xₑ, yₑ) = g₀ / e¹. I know there are multiple solutions to this, but at least one is sufficient.

So far I have

def f(x, y, g0,x0,y0,sigma_x,sigma_y,offset):
    return offset + g0* np.exp(-(((x-x0)**(2)/(2*sigma_x**(2))) + ((y-y0)**(2)/(2*sigma_y**(2)))))

All variables taken as parameters are known as they were extracted from a curve fit.

I understand that taking the derivative in x and setting f() = 0 and similarly in y, gives a solvable linear system for (x,y), but this seems like overkill to manually implement, there must be some library or tool out there that can do what I am trying to achieve?

Answer 1

There are an infinite number of possibilities (or possibly 1 trivial or none in special cases regarding the value of g0). A solution can be computed analytically in constant time using a direct method. No need for approximations or iterative methods to find roots of a given function. It is just pure maths.

Gaussian kernel have interesting symmetries. One of them is the invariance to the rotation when the peak is translated to (0,0). Another on is that the 1D section of a 2D gaussian surface is a gaussian curve.

Lets ignore offset for a moment: it does not really change the problem (it is just a Z-axis translation) and add additional useless term for the resolution.

The geometric solution to is problem is an ellipse so the solution (xe, ye) follows the conic expression : (xe-x0)² / a² + (ye-y0)² / b² = 1. If sigma_x = sigma_y, then the solution is simpler : this is a circle with the expression (xe-x0)² + (ye-y0)² = r. Note that a, b and r are dependant of the searched value and the kernel parameters (eg. sigma_x). Changing sigma_x and sigma_y is like stretching the space, and so the solution similarly. Changing x0 and y0 is like translating the space and so the solution too.

In fact, we could solve the problem for the simpler case where x0=0, y0=0, sigma_x=1 and sigma_y=1. Then we can apply a translation, followed by a linear transformation using a transformation matrix. A basic multiplication 4x4 matrix can do that. Solving the simpler case is much easier since there are are less parameter to consider. Actually, g0 and offset can also be partially discarded of f since it is on both side of the expression and one just need to solve the linear equation offset + g0 * h(xe,ye) = g0 / e so h(x,y) = 1 / e - offset / g0 where h(xe, ye) = exp(-(xe² + ye²)/2). Assuming we forget the translation and linear transformation for a moment, the problem can be solve quite easily:

h(xe, ye) = 1 / e - offset / g0
exp(-(xe² + ye²)/2) = 1 / e - offset / g0
-(xe² + ye²)/2 = ln(1 / e - offset / g0)
xe² + ye² = -2 * ln(1 / e - offset / g0)

That's it! We got our circle expression where the radius r is -2*ln(1 / e - offset / g0)! Note that ln in the expression is basically the natural logarithm.

Now we could try to find the 4x4 matrix coefficients, or actually try to directly solve the full expression which is finally not so difficult.

offset + g0 * exp(-((x-x0)²/(2*sigma_x²) + (y-y0)²/(2*sigma_y²))) = g0 / e
exp(-((x-x0)²/(2*sigma_x²) + (y-y0)²/(2*sigma_y²))) = 1 / e - offset / g0
-((x-x0)²/(2*sigma_x²) + (y-y0)²/(2*sigma_y²)) = ln(1 / e - offset / g0)
((x-x0)²/sigma_x² + (y-y0)²/sigma_y²)/2 = -ln(1 / e - offset / g0)
(x-x0)²/sigma_x² + (y-y0)²/sigma_y² = -2 * ln(1 / e - offset / g0)

That's it! We got you conic expression where r = -2 * ln(1 / e - offset / g0) is a constant, a = sigma_x and b = sigma_y are the unknown parameter in the above expression. It can be normalized using a = sigma_x/sqrt(r) and b = sigma_y/sqrt(r) so the right hand side is 1 fitting exactly with the above expression but this is just some math details.

You can find one point of the ellipse easily since you know the centre of the ellipse (x0, y0) and there is at least 1 point at the intersection of the line y=y0 and the above conic expression. Lets find it:

(x-x0)²/sigma_x² + (y0-y0)²/sigma_y² = -2 * ln(1 / e - offset / g0)
(x-x0)²/sigma_x² = -2 * ln(1 / e - offset / g0)
(x-x0)² = -2 * ln(1 / e - offset / g0) * sigma_x²
x = sqrt(-2 * ln(1 / e - offset / g0) * sigma_x²) + x0

Note there are two solutions (-sqrt(...) + x0) but you only need one of them. I hope I did not make any mistake in the computation (at least the details should be enough to find it easily) and the solution is not a complex number in your case. The benefit of this solution is that it is very very fast to compute.

The final solution is:
(xe, ye) = (sqrt(-2*ln(1/e-offset/g0)*sigma_x²)+x0, y0)