2

我正在使用 python、scipy、numpy 等。

我想根据像素到图像质心的距离来计算灰度图像强度值的直方图。以下解决方案有效，但速度很慢：

import matplotlib.pyplot as plt
from scipy import ndimage
import numpy as np
import math

# img is a 2-dimensionsl numpy array
img = np.random.rand(300, 300)
# center of mass of the pixels is easy to get
centerOfMass = np.array(list(ndimage.measurements.center_of_mass(img)))

# declare histogram buckets
histogram = np.zeros(100)

# declare histogram range, which is half the diagonal length of the image, enough in this case.
maxDist = len(img)/math.sqrt(2.0)

# size of the bucket might be less than the width of a pixel, which is fine.
bucketSize = maxDist/len(histogram)

# fill the histogram buckets
for i in range(len(img)):
    for j in range(len(img[i])):
        dist = np.linalg.norm(centerOfMass - np.array([i,j]))
        if(dist/bucketSize < len(histogram)):
            histogram[int(dist/bucketSize)] += img[i, j]

# plot the img array
plt.subplot(121)
imgplot = plt.imshow(img)
imgplot.set_cmap('hot')
plt.colorbar()
plt.draw()

# plot the histogram
plt.subplot(122)
plt.plot(histogram)
plt.draw()

plt.show()

正如我之前所说，这可行，但速度很慢，因为您不应该在 numpy 中以这种方式对数组进行双循环。有没有更有效的方法来做同样的事情？我假设我需要对所有数组元素应用一些函数，但我也需要索引坐标。我怎样才能做到这一点？目前，一张 1kx1k 的图像需要几秒钟，这非常慢。

4

1 回答 1

2

所有 numpy 分箱函数 ( bincount, histogram, histogram2d... 都有一个weights关键字参数，你可以用它来做一些非常奇怪的事情，比如你的。这就是我的做法：

rows, cols = 300, 300
img = np.random.rand(rows, cols)

# calculate center of mass position
row_com = np.sum(np.arange(rows)[:, None] * img) / np.sum(img)
col_com = np.sum(np.arange(cols) * img) / np.sum(img)

# create array of distances to center of mass
dist = np.sqrt(((np.arange(rows) - row_com)**2)[:, None] +
               (np.arange(cols) - col_com)**2)

# build histogram, with intensities as weights
bins = 100
hist, edges = np.histogram(dist, bins=bins, weights=img)

# to reproduce your exact results, you must specify the bin edges
bins = np.linspace(0, len(img)/math.sqrt(2.0), 101)
hist2, edges2 = np.histogram(dist, bins=bins, weights=img)

两种方法都没有计时，但从从终端运行两种方法的延迟来看，这明显更快。

于 2013-06-21T16:22:34.747 回答