我假设您有一个显示某种峰的图像，并且您有兴趣获取该峰在X和Y方向上的偏度和峰度（以及标准差和质心）。
令人好奇的是，我并没有发现任何python图像分析程序包实现了这些功能。OpenCV有一个moments函数，我们应该能够从这些函数中获得偏度，但moments只到3阶，而我们需要4阶才能得到峰度。
为了使事情更轻松快捷，我认为在X和Y方向上获取图像的投影并找出这些投影的统计信息与使用完整图像找出统计信息在数学上是等价的。在下面的代码中，我同时使用这两种方法并表明对于这个平滑的例子来说它们是一样的。使用一个真实带有噪声的图像，我发现这两种方法也提供了相同的结果，但前提是您手动将图像数据强制转换为float64（导入为float 32，“数字内容”导致结果略有不同）。
以下是一个示例。您应该能够将image_statistics()函数剪切并粘贴到您自己的代码中。希望它对某人有用！
import numpy as np
import matplotlib.pyplot as plt
import time
plt.figure(figsize=(10,10))
ax1 = plt.subplot(221)
ax2 = plt.subplot(222)
ax4 = plt.subplot(224)
#Make some sample data as a sum of two elliptical gaussians:
x = range(200)
y = range(200)
X,Y = np.meshgrid(x,y)
def twoD_gaussian(X,Y,A=1,xo=100,yo=100,sx=20,sy=10):
return A*np.exp(-(X-xo)**2/(2.*sx**2)-(Y-yo)**2/(2.*sy**2))
Z = twoD_gaussian(X,Y) + twoD_gaussian(X,Y,A=0.4,yo=75)
ax2.imshow(Z) #plot it
#calculate projections along the x and y axes for the plots
yp = np.sum(Z,axis=1)
xp = np.sum(Z,axis=0)
ax1.plot(yp,np.linspace(0,len(yp),len(yp)))
ax4.plot(np.linspace(0,len(xp),len(xp)),xp)
#Here is the business:
def image_statistics(Z):
#Input: Z, a 2D array, hopefully containing some sort of peak
#Output: cx,cy,sx,sy,skx,sky,kx,ky
#cx and cy are the coordinates of the centroid
#sx and sy are the stardard deviation in the x and y directions
#skx and sky are the skewness in the x and y directions
#kx and ky are the Kurtosis in the x and y directions
#Note: this is not the excess kurtosis. For a normal distribution
#you expect the kurtosis will be 3.0. Just subtract 3 to get the
#excess kurtosis.
import numpy as np
h,w = np.shape(Z)
x = range(w)
y = range(h)
#calculate projections along the x and y axes
yp = np.sum(Z,axis=1)
xp = np.sum(Z,axis=0)
#centroid
cx = np.sum(x*xp)/np.sum(xp)
cy = np.sum(y*yp)/np.sum(yp)
#standard deviation
x2 = (x-cx)**2
y2 = (y-cy)**2
sx = np.sqrt( np.sum(x2*xp)/np.sum(xp) )
sy = np.sqrt( np.sum(y2*yp)/np.sum(yp) )
#skewness
x3 = (x-cx)**3
y3 = (y-cy)**3
skx = np.sum(xp*x3)/(np.sum(xp) * sx**3)
sky = np.sum(yp*y3)/(np.sum(yp) * sy**3)
#Kurtosis
x4 = (x-cx)**4
y4 = (y-cy)**4
kx = np.sum(xp*x4)/(np.sum(xp) * sx**4)
ky = np.sum(yp*y4)/(np.sum(yp) * sy**4)
return cx,cy,sx,sy,skx,sky,kx,ky
#We can check that the result is the same if we use the full 2D data array
def image_statistics_2D(Z):
h,w = np.shape(Z)
x = range(w)
y = range(h)
X,Y = np.meshgrid(x,y)
#Centroid (mean)
cx = np.sum(Z*X)/np.sum(Z)
cy = np.sum(Z*Y)/np.sum(Z)
###Standard deviation
x2 = (range(w) - cx)**2
y2 = (range(h) - cy)**2
X2,Y2 = np.meshgrid(x2,y2)
#Find the variance
vx = np.sum(Z*X2)/np.sum(Z)
vy = np.sum(Z*Y2)/np.sum(Z)
#SD is the sqrt of the variance
sx,sy = np.sqrt(vx),np.sqrt(vy)
###Skewness
x3 = (range(w) - cx)**3
y3 = (range(h) - cy)**3
X3,Y3 = np.meshgrid(x3,y3)
#Find the thid central moment
m3x = np.sum(Z*X3)/np.sum(Z)
m3y = np.sum(Z*Y3)/np.sum(Z)
#Skewness is the third central moment divided by SD cubed
skx = m3x/sx**3
sky = m3y/sy**3
###Kurtosis
x4 = (range(w) - cx)**4
y4 = (range(h) - cy)**4
X4,Y4 = np.meshgrid(x4,y4)
#Find the fourth central moment
m4x = np.sum(Z*X4)/np.sum(Z)
m4y = np.sum(Z*Y4)/np.sum(Z)
#Kurtosis is the fourth central moment divided by SD to the fourth power
kx = m4x/sx**4
ky = m4y/sy**4
return cx,cy,sx,sy,skx,sky,kx,ky
#Calculate the image statistics using the projection method
stats_pr = image_statistics(Z)
#Confirm that they are the same by using a 2D calculation
stats_2d = image_statistics_2D(Z)
names = ('Centroid x','Centroid y','StdDev x','StdDev y','Skewness x','Skewness y','Kurtosis x','Kurtosis y')
print 'Statistis\t1D\t2D'
for name,i1,i2 in zip(names, stats_2d,stats_pr):
print '%s \t%.2f \t%.2f'%(name, i1,i2)
plt.show()
输出结果：
Statistis 1D 2D
Centroid x 99.48 99.48
Centroid y 99.86 99.86
StdDev x 19.83 19.83
StdDev y 9.92 9.92
Skewness x 0.00 0.00
Skewness y 0.00 0.00
Kurtosis x 3.00 3.00
Kurtosis y 3.00 3.00
同样，我还提供了一个2D版本的函数image_statistics_2D()，它应该与1D版本提供相同的结果。
最后，我还要提一下，根据您对图像的确切操作，您可能会考虑使用ImageJ进行图像分析——但是要当心！Moments插件可以让您计算偏度、峰度等。ImageJ在Analyze>>Set Measurements菜单中确实有“skewness”和“kurtosis”，但我认为这
