大约12年前，我实施了一个算法roifill在图像处理工具箱中。在该算法中，将替换(或填充)一组指定的像素，使每个替换的像素等于其北、东、南和西邻居的平均值。我执行了一个迭代算法，虽然有效，但速度很慢。当我问MATLAB的创造者和公司的联合创始人克里夫如果他知道更好的方法，他立即回答说:“对我来说，它看起来像一个稀疏线性系统。为什么不直接用反斜杠呢?”

头一巴掌!当然可以。

现在，我又看到了一组非常相似的等式，并再次考虑使用稀疏反斜杠。

在Tarboton的论文中，上坡面积的计算是递归的。也就是说，一个像素的上坡面积等于1加上它相邻像素上坡面积的一部分。(比分数是相邻水流方向的函数。)这是一个线性方程组。让我们看看如何在MATLAB中建立和求解这个系统。

我们将从创建一个测试“DEM”开始。我喜欢用山峰．

E =山峰;imshow (E, [],“InitialMagnification”，“健康”）

接下来，使用函数计算像素流方向dem_flow从我的上坡区域贡献的MATLAB中央文件交换．

R = dem_flow(峰);

稀疏系统的每一行和每一列都对应于一个单独的像素。它将方便地使一个矩阵，同样大小的E，“数字”或标签的个别像素的E，这些标签将被用作稀疏矩阵的行和列坐标。

pixel_labels = 1:元素个数(E);pixel_labels =重塑(pixel_labels，大小(E));pixel_labels(1:4,书1:6)

Ans = 150 99 148 197 246 2 51 100 149 198 247 3 52 101 150 199 248 4 53 102 151 200 249

让我们考虑一下北方的邻居。这是一组像素与北部邻居:

Pixels_with_north_neighbors = pixel_labels(2:结束，:);

这里是对应的北方邻居

North_neighbors = pixel_labels(1:end-1，:);

从像素的北邻向内的角度是- /2:

inward_angle = -π/ 2;

一个像素的北邻居有一个指向特定方向的流。这个方向包含在R．如果气流方向与内角相同，那么来自邻居的气流比例为1.0。如果气流方向指向距离内角大于/4弧度的方向，那么来自邻居的气流的比例是0.0。对于之间的角度差异，分数从1.0到0.0线性变化。

我已经写了简单的函数来计算角差和相应的权重。角差的计算是这样的:

D = mod(1 - 2 + 2) -;

方向权值的计算是这样的:

W = max(1 - abs(d)*4/pi, 0);

通过这些函数，我们可以一次计算所有像素的北邻权值:

w = directional_weight (angular_difference (R (north_neighbors),.．.inward_angle));

暂时忽略其他邻居，上坡面积方程是这样的:

Ap = 1 + wAn

其中，Ap为像素的上坡区域，An为北邻的上坡区域，w为根据北邻的流向确定的方向权值。这样重写方程，看看如何构造线性方程矩阵:

Ap - wAn = 1

所以矩阵的每一行在对角线上都有一个1，它将有负的方向流动权值对于它的每一个邻居。如果我们只使用上面计算的北方权值，我们可以这样构造稀疏矩阵:

i = pixels_with_north_neighbors (:);j = north_neighbors (:);w = w (:);添加对角项Ii = [i;pixel_labels (:));Jj = [j];pixel_labels (:));[w =];(元素个数(pixel_labels), 1)];并构造线性系统矩阵，调用稀疏T =稀疏(ii, jj, ww);间谍(T)

我们必须放大才能清楚地看到非对角线项:

[0 20] ([0 20])

我要打电话了T的流矩阵．完成计算T，我们必须为其他七个相邻方向重复上面的一些步骤，在这里重复太单调了。我就用flow_matrix函数在我MATLAB中央文件交换提交：

T = flow_matrix(E, R);

现在我们可以通过求解矩阵方程来计算每个像素的上坡面积:

Ta = b

其中T是我们刚刚形成的稀疏矩阵，a是每个像素的上坡面积的向量，b是1的向量。

a = T \ ones(numel(E)， 1);A = shape(A, size(E));imshow ([],“InitialMagnification”，“健康”）

我发现使用log(A)可以很好地可视化上坡面积矩阵:

imshow(日志(A)、[]“InitialMagnification”，“健康”）

较亮的值有最高的上坡面积。

在本系列的下一篇文章中，我将重新讨论处理平台的问题。

获取MATLAB代码

与MATLAB®7.4发布