gpufit / Gpufit

GPU-accelerated Levenberg-Marquardt curve fitting in CUDA
MIT License
312 stars 92 forks source link

Use GAUSS_2D_ROTATED to fit a 2D Gaussian, the returned state is 2 #98

Closed PengfeiRen96 closed 3 years ago

PengfeiRen96 commented 3 years ago

The code is: import numpy as np import matplotlib.pyplot as plt import pygpufit.gpufit as gf

def gauss2d_angle(x, y, amp, x0, y0, sigma_x, sigma_y, theta): a = (np.cos(theta) 2) / (2 sigma_x 2) + (np.sin(theta) 2) / (2 sigma_y 2) b = -(np.sin(2 theta)) / (4 sigma_x 2) + (np.sin(2 theta)) / (4 sigma_y 2) c = (np.sin(theta) 2) / (2 sigma_x 2) + (np.cos(theta) 2) / (2 sigma_y 2) out = amp np.exp(- (a ((x - x0) * 2) + 2 b (x - x0) (y - y0)

if name == 'main':

# cuda available checks
print('CUDA available: {}'.format(gf.cuda_available()))
if not gf.cuda_available():
    raise RuntimeError(gf.get_last_error())
print('CUDA versions runtime: {}, driver: {}'.format(*gf.get_cuda_version()))

# number of fit points
size_x = 32
number_points = size_x * size_x

# average noise level
average_noise_level = 0.01

# initialize random number generator
np.random.seed(0)

# tolerance
tolerance = 0.0001

# maximum number of iterations
max_number_iterations = 100

# generate x and y values
g = np.arange(size_x)
yi, xi = np.meshgrid(g, g, indexing='ij')
xi = xi.astype(np.float32)
yi = yi.astype(np.float32)
n_fits = 1

# generate data
data = gauss2d_angle(xi, yi, 0.4, 16, 16, 1.4, 2.5, np.pi/4)[np.newaxis, :]

initial_parameter = np.array((0.4, 14, 17, 1, 1, 1, 0.1), dtype=np.float32)[np.newaxis, :]

# run Gpufit
parameters, states, chi_squares, number_iterations, execution_time = gf.fit(data, None, gf.ModelID.GAUSS_2D_ROTATED,
                                                                            initial_parameter, tolerance,
                                                                            max_number_iterations)
data_fitted = gauss2d_angle(xi, yi, *parameters[0, :-2], parameters[0, -1])

fig, ax_list = plt.subplots(1, 2, figsize=(4, 8))
ax_list[0].imshow(data[0].reshape(size_x, size_x), interpolation='nearest')
ax_list[1].imshow(data_fitted.reshape(size_x, size_x), interpolation='nearest')
fig.tight_layout()
plt.show()
jkfindeisen commented 3 years ago

If the return state is really 2 than that is strange because according to the documentation it should only be 0 or -1 indicating that there was an error or not. Maybe you mean the fit state, which can be 2 and means that the Hessian is singular. That might hint that the used fit model is not a good model here, some parameters may not be needed or something else is wrong.

superchromix commented 3 years ago

This could happen if the input data is perfectly symmetric. In this case, the rotation angle has no effect on chi-square.

PengfeiRen96 commented 3 years ago

@superchromix Thanks for your help, the problem was solved when I modified the initial parameter.