2D Gaussian curve fit function

[minhviet0308]

@isaachroberts43 Can you figure out a way to create a function in python that does the 2D Gaussian curve fit for a matrix like the IDL code? Apparently, python does not have a built-in function or method that does this. The function for it in IDL is gauss2dfit().

[isaachroberts43]

@minhviet0308 I talked to Suman about implementing this and we found documentation for the IDL function https://www.l3harrisgeospatial.com/docs/gauss2dfit.html, I think if we implement this function that would work. I can help you tomorrow in person.

[sumanchapai]

@minhviet0308 @isaachroberts43 what's the progress on this?

[minhviet0308]

@sumanchapai I completed the Gaussian 2d curve fit using scipy. There are two things that needs to be done:

1. I need to compare the results with the IDL version to see if they match.
2. It takes a while to go through one iteration of matrices that we're interested in (in this case I suppose it's the hot pixels), so I will need help optimizing this code.

Here is a demo code with a simple 5x5 matrix. I will need to formulate this into a function for us to use.

import numpy as np 
from scipy.optimize import curve_fit 
import math

M = np.ones((5,5))
M[1:3,1:3] = 4
M[2,2]=10
print(M)

data = []
for index, z in np.ndenumerate(M): 
    data.append((index[0], index[1],z))

A = np.array(data)

def func(data, A0, A1, a, b, h, k): 
    x = data[:, 0]
    y = data[:, 1]
    return A0 + A1*np.exp(-(((x-h)/a)**2+((y-k)/b)**2)/2)

guess = (1,9,2,2,2,2)
params, pcov = curve_fit(func, A[:,:2], A[:,2], guess, maxfev = 5000)

A0 = params[0]
A1 = params[1]
a = params[2]
b = params[3]
h = params[4]
k = params[5]

for index, z in np.ndenumerate(M): 
    x = index[0]
    y = index[1]
    new_z = A0 + A1*np.exp(-(((x-h)/a)**2+((y-k)/b)**2)/2)
    M[x,y] = new_z

print(M)