使用GDAL和Python的最小距离算法

Question

我正在尝试使用GDAL和Python实现最小距离算法进行图像分类。在计算样本区域的平均像素值并将它们存储到一个数组列表（“sample_array”）之后，我将图像读入名为“values”的数组中。使用以下代码循环遍历该数组：

values = valBD.ReadAsArray()

# loop through pixel columns
for X in range(0,XSize):

    # loop thorugh pixel lines
    for Y in range (0, YSize):

        # initialize variables
        minDist = 9999
        # get minimum distance
        for iSample in range (0, sample_count):
            # dist = calc_distance(values[jPixel, iPixel], sample_array[iSample])

            # computing minimum distance
            iPixelVal = values[Y, X]
            mean = sample_array[iSample]
            dist = math.sqrt((iPixelVal - mean) * (iPixelVal - mean)) # only for testing

            if dist < minDist:
                minDist = dist
                values[Y, X] = iSample

classBD.WriteArray(values, xoff=0, yoff=0)

这个过程对于大图片来说非常耗时。因此，我想问一下是否有人知道更快的方法。我不太了解Python中不同变量的访问速度。或者也许有人知道我可以使用的库。 提前感谢， Mario

- Mario Härtwig

---

在纯Python中逐像素处理任何内容都可能会很慢。你可以查看Python Imaging Library（PIL）的功能是否有用。 - Thomas K

Answer 1

你一定要使用NumPy。我处理的一些相当大的栅格数据集，NumPy可以轻松应对。在我的机器上，对于下面的代码，1000 x 1000数组没有任何明显的延迟。下面是这个代码如何工作的解释。

import numpy as np
from scipy.spatial.distance import cdist

# some starter data
dim = (1000,1000)
values = np.random.randint(0, 10, dim)

# cdist will want 'samples' as a 2-d array
samples = np.array([1, 2, 3]).reshape(-1, 1)

# this could be a one-liner
# 'values' must have the same number of columns as 'samples'
mins = cdist(values.reshape(-1, 1), samples)
outvalues = mins.argmin(axis=1).reshape(dim)

cdist()函数计算每个元素与samples中的所有元素之间的“距离”。这将生成一个1,000,000 x 3的数组，其中每行n都具有从原始数组中的像素n到每个样本值[1, 2, 3]的距离。使用argmin(axis=1)函数可以获取每行中最小值的索引，这正是您需要的。快速重塑即可获得图像的矩形格式。

- Robin Kraft

Answer 2

同意Thomas K的观点：使用PIL，或者编写一个C函数并使用ctypes进行封装，或者至少使用一些numPy矩阵操作。

否则，在现有代码上使用pypy（JIT编译的代码在图像代码上可以快100倍）。尝试pypy并告诉我们你得到了多少加速。

底线：绝不要在cPython中本地像这样逐像素地处理内容，解释和内存管理开销会让你崩溃。

- smci