2016-04-06 49 views
2

私は暗い画像(生のフォーマット)を持っており、画像と画像の分布をプロットしています。ご覧のように、16時にピークがありますが、無視してください。私はこのヒストグラムを通してガウス曲線に収まるようにしたい。私はこの方法を以下のように適合させました: Un-normalized Gaussian curve on histogramしかしながら;私のGaussian fitは決して近くに来ない。私はプロットのための正しいフォーマットにイメージを変えることで何か間違っているのですか、それとも間違っているのでしょうか? enter image description here Gaussian distribution of the imageヒストグラムpythonの正規化されていないガウス分布を合わせる

これは、私は、このデータを生成するために使用する現在のコードです:

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

def fitGaussian(x,a,mean,sigma): 
    return (a*np.exp(-((x-mean)**2/(2*sigma)))) 

fname = 'filepath.raw' 
im = np.fromfile(fname,np.int16) 
im.resize([3056,4064]) 

plt.figure() 
plt.set_cmap(viridis) 
plt.imshow(im, interpolation='none', vmin=16, vmax=np.percentile(im.ravel(),99)) 
plt.colorbar() 
print 'Saving: ' + fname[:-4] + '.pdf' 
plt.savefig(fname[:-4]+'.pdf') 

plt.figure() 
data = plt.hist(im.ravel(), bins=4096, range=(0,4095)) 

x = [0.5 * (data[1][i] + data[1][i+1]) for i in xrange(len(data[1])-1)] 
y = data[0] 

popt, pcov = curve_fit(fitGaussian, x, y, [500000,80,10]) 
x_fit = py.linspace(x[0], x[-1], 1000) 
y_fit = fitGaussian(x_fit, *popt) 
plt.plot(x_fit, y_fit, lw=4, color="r")   

plt.xlim(0,300) 
plt.ylim(0,1e6) 
plt.show() 

EDIT:(ルブロション仮面劇への応答)

私はまだ16でビンを削除した場合同じフィットを得る: enter image description here

+0

私はあなたが何をしているのかよく分かっていませんが、ヒストグラムのスパイク20を取り除きたい場合があります。 –

+0

@ReblochonMasque 16でビンをなくしても、結果は、右にシフトしただけです。私の編集を参照してください。 – SjonTeflon

+1

私はゴーセージをすべてのビンに合わせるので、このカーブを得ていると思います。大部分はゼロです。あなたはカーブをゼロでないビンにのみフィットさせることができますが、これは面倒です。 – kazemakase

答えて

2

適合したガウスはあまりにも低いように見える大部分がゼロであるすべてのビンに収まるからです。解決策は、ガウス関数を非ゼロビンにのみフィットさせることです。

plt.histの代わりにnp.histogramを使用してビンビンを取得しますが、これは味の問題です。重要な部分はxhyh

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

# Generate example data 
x = np.random.randn(100000) * 50 + 75 
x = np.round(x/10) * 10 
x = x[x >= 20] 

yhist, xhist = np.histogram(x, bins=np.arange(4096)) 

xh = np.where(yhist > 0)[0] 
yh = yhist[xh] 

def gaussian(x, a, mean, sigma): 
    return a * np.exp(-((x - mean)**2/(2 * sigma**2))) 

popt, pcov = curve_fit(gaussian, xh, yh, [10000, 100, 10]) 

plt.plot(yhist) 
i = np.linspace(0, 300, 301) 
plt.plot(i, gaussian(i, *popt)) 
plt.xlim(0, 300) 

enter image description here

P.S.の定義ですシグマは通常、分散ではなく標準偏差を表します。だから私はgaussianの関数でそれを二乗したのです。

関連する問題