Python3, curve_fitによる2次元画像データのガウシアンフィッティング

Question

###前提・実現したいこと
Python3とscipy.optimize.curve_fitを用いて二次元画像データからガウシアンフィッティングを行いたいのですが,

Result from function call is not a proper array of floats.

のエラーメッセージの対処がわからず困っています。
対処法をお教えくださると幸いです。

###該当のソースコード

import numpy as np
from PIL import Image as pil
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

data = np.array(pil.open("figs/image.png"))
data = data/np.max(data)

x,y=np.meshgrid(np.linspace(0,data.shape[1],data.shape[1]),np.linspace(0,data.shape[0],data.shape[0]))

def twoDgaussian(X,wx,wy,x0,y0):
    x,y=X
    z=np.exp(-(x-x0)**2/wx**2-(y-y0)**2/wy**2)
    return z

initial=(50,50,500,500)

popt,pcov=curve_fit(twoDgaussian,(x,y),data,initial)

###発生している問題・エラーメッセージ

-------------------------------------------------------------------------
ValueError                              Traceback (most recent call last)
ValueError: object too deep for desired array

-------------------------------------------------------------------------
error                                   Traceback (most recent call last)
<ipython-input-133-8afc316156ee> in <module>()
----> 1 popt,pcov=curve_fit(twoDgaussian,(x,y),data,initial)

C:\Users\***\Anaconda3\lib\site-packages\scipy\optimize\minpack.py in curve_fit(f, xdata, ydata, p0, sigma, absolute_sigma, check_finite, bounds, method, jac, **kwargs)
    734         # Remove full_output from kwargs, otherwise we're passing it in twice.
    735         return_full = kwargs.pop('full_output', False)
--> 736         res = leastsq(func, p0, Dfun=jac, full_output=1, **kwargs)
    737         popt, pcov, infodict, errmsg, ier = res
    738         cost = np.sum(infodict['fvec'] ** 2)

C:\Users\***\Anaconda3\lib\site-packages\scipy\optimize\minpack.py in leastsq(func, x0, args, Dfun, full_output, col_deriv, ftol, xtol, gtol, maxfev, epsfcn, factor, diag)
    385             maxfev = 200*(n + 1)
    386         retval = _minpack._lmdif(func, x0, args, full_output, ftol, xtol,
--> 387                                  gtol, maxfev, epsfcn, factor, diag)
    388     else:
    389         if col_deriv:

error: Result from function call is not a proper array of floats.

###追加点

def twoDgaussian(X,wx,wy,x0,y0):
    x,y=X
    z=np.exp(-(x-x0)**2/wx**2-(y-y0)**2/wy**2)
    return z.ravel()

initial=(50,50,500,500)
data_ravel=data.ravel()
popt,pcov=curve_fit(twoDgaussian,(x,y),data_ravel,initial)

とすることで通るようになったものの、出力結果が
popt=array([ 1.05774768e+00, 1.23219802e-02, 5.00141975e+02, 5.00041068e+02])
でwxとwyは同程度のはずなのでまだおかしい。

Comments

[退会済みユーザー]

どの変数が適切なfloatになっていないか呈示できそうでしょうか？何行目の何がどうこう、というエラーメッセージが表示されていなさそうという理由です。

[oknd1]

すみません。表示されているメッセージとしては他にないのではっきりとわからないのですが、関数twoDgaussianの出力のことを指していると思います。ただ、この関数の出力は確認するとfloatのarrayになっているようでした。

[退会済みユーザー]

Q.1 curve_fit(以下略）で twoDgaussian(X,wx,wy,x0,y0)　を呼び出す際に、twoDgaussian()の引数が省略されているのはこれで大丈夫でしょうか？
Q.2 差し支えなければimage.pngをアップロードしてください

[oknd1]

アップロードしました。引数の省略ですが、チュートリアルなどでも省略しているのと、一変数関数のフィッティングを試してみて問題なく通ったことからそこは問題ないと思います。第一引数を独立変数として扱ってそれ以外はパラメータとして扱われるようです。

Answer 1

完成形がよく分からないので、回答でありながら確認になってしまいますが、
以下のコード添付の画像で動くようにするイメージでしょうか？

Stackoverflowのコードを動くようにしただけですが...

import scipy.optimize as opt
import numpy as np
import pylab as plt
from PIL import Image as pil

#define model function and pass independant variables x and y as a list
def twoD_Gaussian(XY, amplitude, xo, yo, sigma_x, sigma_y, theta, offset):
    x,y = XY[0:2]
    xo = float(xo)
    yo = float(yo)
    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)
    g = offset + amplitude*np.exp( - (a*((x-xo)**2) + 2*b*(x-xo)*(y-yo)
                            + c*((y-yo)**2)))
    return g.ravel()


# Create x and y indices
x = np.linspace(0, 200, 201)
y = np.linspace(0, 200, 201)
x, y = np.meshgrid(x, y)

#create data
data = twoD_Gaussian((x, y), 3, 100, 100, 20, 40, 0, 10)


# plot twoD_Gaussian data generated above
plt.figure()
plt.imshow(data.reshape(201, 201))
plt.colorbar()

# add some noise to the data and try to fit the data generated beforehand
initial_guess = (3,100,100,20,40,0,10)

data_noisy = data + 0.2*np.random.normal(size=data.shape)

popt, pcov = opt.curve_fit(twoD_Gaussian, (x, y), data_noisy, p0=initial_guess)

data_fitted = twoD_Gaussian((x, y), *popt)

fig, ax = plt.subplots(1, 1)
ax.hold(True)
ax.imshow(data_noisy.reshape(201, 201), cmap=plt.cm.jet, origin='bottom',
    extent=(x.min(), x.max(), y.min(), y.max()))
ax.contour(x, y, data_fitted.reshape(201, 201), 8, colors='w')
plt.show()

Comments

[oknd1]

slash様

まさにそうです。ただstackoverflowの提示していただいたページは参照にしていたもののpython2.7系とpython3系で微妙に違いがあるようで苦戦していました。
先程ravel()を適宜加えることで通るようにはなったのですが、まだ出力結果がおかしい状態です。
ありがとうございます。

[oknd1]

出力結果がおかしいといった話でしたが、別にテストデータを生成して確認してみたところ問題ないようでしたのでフィッティングの流れとしてはこれで問題ないようです。
ほぼ解決しましたありがとうございました。