使用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

Answer 1

我同意Thomas K的看法：可以使用PIL库，或者自己写一个C语言的函数，然后用比如ctypes来包装它。至少也可以使用一些numPy的矩阵操作。

另外，你也可以在现有代码上使用pypy（JIT编译的代码在处理图像时可以快100倍）。试试pypy，然后告诉我们你得到了多大的速度提升。

总之：在cPython中不要像这样逐个像素地处理图像，解释和内存管理的开销会让你受不了。

Answer 2

你绝对应该使用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()这个函数用来计算values中每个元素到samples中每个元素的“距离”。这样就生成了一个1,000,000 x 3的数组，其中每一行n表示原始数组中第n个像素到每个样本值[1, 2, 3]的距离。argmin(axis=1)会给你每一行中最小值的索引，这正是你想要的。然后快速调整一下形状，就能得到你期望的图像的矩形格式。