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

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?

CodePudding user response：

You can use the scipy.optimize.fsolve function to find the (xₑ, yₑ) values that satisfy the equation f(xₑ, yₑ) = g₀ / e¹. This function uses a numerical root-finding algorithm to find the roots of a system of non-linear equations.

Here's an example of how you can use scipy.optimize.fsolve to find the (xₑ, yₑ) values that satisfy the equation f(xₑ, yₑ) = g₀ / e¹:

from scipy.optimize import fsolve

# Define the function f(x, y)
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)))))

# Define the function that we want to find the roots of
def g(xy, g0, x0, y0, sigma_x, sigma_y, offset):
    x, y = xy
    return f(x, y, g0, x0, y0, sigma_x, sigma_y, offset) - g0 / np.exp(1)

# Define the initial guess for the root
x0_guess = x0
y0_guess = y0

# Find the roots of the function g(x, y)
x, y = fsolve(g, (x0_guess, y0_guess), args=(g0, x0, y0, sigma_x, sigma_y, offset))

# Print the result
print(f"x = {x}, y = {y}")

In this example, fsolve will use the initial guess (x0_guess, y0_guess) as a starting point and iteratively try to find the roots of the function g(x, y). If the function g(x, y) has multiple roots, fsolve will return only one of them (the one closest to the initial guess).

CodePudding user response：

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)